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I. INTRODUCTION 


A. BACKGROUND 

The NPS Operations Research Department has made Radio Shack TRS-80 
Model 100 computers available to a limited number of students. Experience has shown 
that these students have made relatively limited use of these computers. The main 
easons seem to be: 

e This computer uses the BASIC programming language. NPS OA — do 
noi receive instruction in this language but_aré required to learn FORTRAN and 
APL. Although FORTRAN and BASIC have’ similar structures. most OA 
students believe they do not have tme to learn a third computer language. 


e Only a few M-100 programs are currently available at NPS that directly relate to 
course work. 


e Whiutng — debugging programs for the M100 can take a considerable amount 
of umé. Many OA Sttide eis believe that time would be more usefully spent 
pursuing other approaches to their studies. 

° The M100 has thus far not been distributed to all OA “7 Therefore. 

fessors have not been able to require students to use the M100. It has been 


icieaced to “nice to have” status instead of being included as an essential 
teacning tool. 


B. GENERAL | 

The purpose of this paper is to develop a collection of programs in BASIC for 
the M100 that students can use in the combat modeling courses of the OA curriculum. 
The programs make extensive use of subroutines which allow students to run programs 
“off the shelf” or build their own programs from the subroutines. The programs are 
extensively documented so that students who learn BASIC can use the printed 
programs as study aids to understand the algorithms involved. Some of the programs 
are tutorial. 

In testing situations professors are often limited to problems that are 
arnithmetically very simple. If programs for the M100 were available, professors would 
be able to give more intricate test problems without placing undue emphasis on manual 
arithmetic calculation. 

Additionally, when students leave NPS, they will be able to take with them a 


series of programs with which they have grown familiar during their course of study. 


I. FEATURES COMMON TO MORE THAN ONE PROGRAMI 


A. GENERAL CHARACTERISTICS OF THE MODEL 100 

The M100 ts a versatile and very portable computer. As provided to Naval 
Postgraduate School students, it has either 24K or 32K of RAM for storage of 
variables during program execution and for storage of programs and other files. 
Programs and files remain active while the computer is turned off. There is an internal 
300 baud modem which facilitates transfering programs to other computers for storage. 
The eight line LCD screen limits the graphics display capability of the M100. 
However, output in character form can be written to RAM files and reviewed after 
program execution is complete. The version of BASIC used in the M100, creation of 
text files, use of the modem, etc. are explained in Reference 1. This thesis assumes that 
readers are somewhat familiar with Reference 1. 

Although the programs have statements which print intermediate and final results 
to the screen, the operator may wish to check the status of a variable that is not 
printed by the program. To do this | 

e Stop the program by depressing the SHIFT and BREAK keys simultaneously. 


° wipe a BASIC statement to print the desired variables: and hit the ENTER 
utton. 


e Start the program again at the place it stopped by typing CONT and hitting the 
ENTER button. P pped by typing CO! g 


Most of the programs in this thesis do not include statements to end the 
program. This is because the programs cycle back to the start of the program allowing 
the operator to enter a new set of parameters without restarting the program. To end 
the program 

e Depress the SHIFT key and the BREAK key at the same time. 
e To rerun the program from the start, depress the F4 key. 


e To return to the main menu, depress the F8 key. 


B. COMMON TERMINOLOGY 
BASIC variables are refered to in the text using capital letters. Since M100 
BASIC only differentiates between variables based on the first two letters of the 


variable name, most BASIC variables in the following programs are combinations of 


two capital letters, e.g. BA, CR, FT. Names of matrices are also specified using capital 
letters. BASIC permits matrices to be dimensioned using variables. If, for instance, 
AB=3 and AC=4, then the BASIC statement DIM ZZ(AB,AC) would dimension a 
three by four matrix named ZZ. When capital letters are used inside the parentheses, 
the size of the matrix is being specified. When a small letter is used inside the 
parentheses, a particular element of the matrix is being specified. For example, ZZ(i.}) 
refers to the element of ZZ in the ith row and jth column. When a matrix has more 
than one dimension, the first number from the right is refered to as the column number 
and the second number from the right as the row number. 

When the mathematical theory behind a program is discussed, the variables used 


will be small letters or Greek letters. 


C. ENPOT FILES 

All of the programs presenied in this thesis permit data to be entered from a text 
file. These text files must be created using the M100’s text editor before the program is 
run. The name of the text file for each program is similar to the name of the program 
It supports and is given in the documentation for each program under the section 
titled, “Input”. The contents of the input file are also explained in the appropnate 
“Input” section. 

Numbers in input files must be separated by a space, comma, or return. A 
comma and a return may not be placed together without a number between them. 
Otherwise, the M100 will enter an extra zero at that point. A comma and a return— 
together also cause data elements which follow to be in the wrong places in their 
matrices and/or wrong numbers to be read as matnx dimensions. Do not put blank 


lines’ in the data for the same reason. 


D. FORMULA TOKENIZATION 

There are some programs which must be adapted to use different equaiions at 
the same point in the program depending on the application. For example, the 
detection theory simulation in Chapter 3 must be able to handle many functions for 
the probability of detection of a sensor. When the function needs to be changed, the 


operator may: 


! Two returns together. 
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e Stop the program, call the BASIC editor for the lines that need to be changed, 
and restart the program after the function has been edited, or 


e Use a tokenization routine to change the function without stopping the program. 


This section describes a subroutine (see Figure 2.1) which performs that tokenization. 


'ToKenize DF 
BS="DF="4+DFS+CHRS$(0 } 


'Tokenize/execute B$ 
BO=VARPTR(BS ):B1=PEEK(BO+1 )+256*PEEK( BO+2 ):CALL1606,0,Bl1 
CALL2499,0,63105:RETURN 





Figure 2.1 Formula Tokenization Section. 


That subroutine is taken from [Ref. 2:pp. 58-60]. Tokenization converts a string 
variable into an executable BASIC statement. In the example in Figure 2.1, the right 
hand side of the function equation was stored as a string variable, DFS before this 
subroutine was called. For example, if the equation was DF=X+Y, then DFS is the 
string, “X+Y”. Line 1410 adds an equal sign and appropriate left hand side to DFS to 
form BS, a complete assignment statement in string format. In the example above line 
1410 adds *DF=” to DFS to form BS. Line 1451 computes memory location, Bl, of 
BS and calls the machine level subroutine in the M100 read only memory (ROM) 
beginning at memory location 1606. Subroutine 1606 converts the BASIC keywords in 
BS into their one byte tokens and then stores them in executable form at memory 
location 63105. Line 1455 calls the machine level subroutine beginning at memory 


location 2499 which executes the statement stored at memory location 63105. 
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It]. DETECTION SIMULATION PROGRAM 


A. GENERAL 

The problem addressed by this program is estimating the probability of detecting 
a target whose location is given by a bivariate normal distribution when the location of 
each detecting sensor also has a separate BVN distribution. The target distribution is 
BVN(0,0,07,,07y,Px y)- The location of each sensor, S,, is distributed BVN(H,, , Hy.» 
oy. o v Pu.v All distribution parameters are stored in an input text file. The 
program models two sensor types: 

e A Deterministic Sensor. This sensor has detection probabilities of 1, 0 and 1 in 
three concentric circular bands around the sensor's location. An example of this 
type of sensor would be a sonobouy which detects all targets out to range rj, 
detects no targets between ranges r, and rp», detects all targets in a convergence 
zone between r4 and r3, and detects no targets beyondr3 whereO Sry Sry S 
r3. 

e A Probabilistic Sensor. This sensor has a continuous, radially symmetric 
detection pattern. The probability of detection, Py, is a function of range and 
may also be a function of one or more other parameters. The default function, 
D(r,6), for calculating P, is a Carlton function where b is a scaling parameter 
(See Equation 3.1). Figure 3.1 shows graphs of several Carlton detection 


functions. 


Daestd 
Poem) = <7! (22 (eqn 3.1) 


For either type of sensor the program calculates the overall probability of detection 
using D(r). The deterministic sensor portion of the program can also be used to 
calculate a “cookie cutter” approximation of the actual detection function. The 
cookie-cutter approximation has the form D(r)=1 when rSr, and D(r)=0 when 


r>r,, for some specified detection range, r,. The radius, r, of the cookie-cutter 


* 
approximation is the lethal radius of the actual detection function for a sensor. Lethal 
radius is a term borrowed from the damage functions of firing theory. It is used here 


to help readers who are familiar with firing theory, not because being detected by a 
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sensor is inherently lethal. The lethal radius is a scalar measure of a sensor's detection 


ability and is computed using Equation 3.2. See {Ref. 3:pp. 12-15]. 
Lethal Radius = C2(> rD(r) dr]? (eqimse: ) 


To calculate a cookie-cutter approximation, the lethal radius of D(r) must be computed 
off line and included as a sensor parameter.* For example, the lethal radius of the 


default Carlton detection function 1s 2b, 


« THREE CARLTON DETECTION FUNCTIONS 
with b= 1, 3, = 
© 


to | 


ae 
G 
= 
Ole) 
; © 
= 
ob) 
ox 





Figure 3.1 Graph of Carlton Detection Functions. 


Closed form solutions are available for some specific cases that can also be 
computed with this program, e.g. a group of sensors with Carlton detection functions 
at fixed locations. The program will verify these closed form solutions. However, it 
will also estimate probability of detection for problems without closed form solutions, 
e.g. a group of cookie-cutter sensors in which each sensor has a different BVN location 


distribution. 


h ore the second example in the section of example problems at the end of this 
chapter. 
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The program uses a Monte Carlo simulation to estimate Py. Because the Monte 


Carlo simulation is time consuming, faster numerical integration approximations are 


also available for the deterministic sensor. [lowever, the particular numerical integration 


technique used in this program will not produce correct results if sensor areas of coverage 


overlap. More sophisticated numerical integration techniques are available for special 


cases of the probabilistic sensor, but are not included in this program. 


B. 


EXPLANATION OF VARIABLES 

AL is the probability that the computed confidence level does not contain the 
true probability of detection for the associated standard normal CDF value. 
A2(6,6) 1s the intermediate calculation matrix for the Romberg integration 
subroutine. 

BO and BI are the memory location parameters of BS. 

BS is the string that holds the entire equation for DF in the tokenization process. 
D1 and D2 are intermediate calculation values. 

DF is the value of the detection function. 

DFS is the string that hold the right hand side of the equation for DF. 

F is the value of the function being integrated at a specific point. 

Fl is, in the data input section (lines100-150), the selection variable indicating 
whether all sensors have the same parameters. If Fl=1, then sensor parameters, 
other than the five for location, will be listed only once: after the location 
parameters of S,. If Fl=0, then values of all parameters for all sensors must be 
entered explicitly in DSIN.DO. In the rest of the program, beyond line 200, F1 
is a flag determining whether the calculation is based on the deterministic 
(F1l= 1) or the probabilistic (Fl = 2) sensor. 

F2 is the selection variable for whether all sensors locations are distributed with 
= o,=0. 1 - yes; 2 - no. 

F3 is the selection variable for whether a Monte Carlo simulation (F3=1) or a 
numerical integration (F3= 2) is used. 

H is the radius of circular limits of integration. 

I, 11, 12 are loop counters. | 

IN is the value of the integral from a numerical approximation. 

J1 through J9 are loop counters. 

NS is the number of sensors. 


NP is the number of parameters for each sensor in addition to location. 
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NR is the number of repetitions in a Monte Carlo simulation. 

PD is the probability of detection. 

RF is (1-RH2)>. 

RH is the correlation between components of the BVN target location 
distribution. 

Sl and S2 are standard deviations of the components of the target location 
distribution. 

TE is a temporary storage variable. 

TIS and T2$ are times at beginning and end of a calculation. 

Tl is the time elapsed during a calculation. 

Vil and V2 are variances of the components of the target location distribution. 
XS and XT are specific values of the first component, i.e. the mean X position, 
of the BVN distributions of a sensor or the target respectively. 

X(NS,5+ NP) is the parameter matrix for sensors. There is one row per sensor 
and one column per parameter. Columns one through five are for the means, 
Hus. and HyS. the standard deviations, Cus. and OVS. and the correlation, 
Pco., of the location of each sensor, S;. These parameters are in the same 
coordinate system as the target location distribution. 

YS and YT are specific values of the second component, i.e. the mean Y position, 
of the BVN distributions of a sensor or the target respectively. 


Zl and Z9 are selection variables. 


INPUT 
1. Input File 
Before this program 1s run, a text file, DSIN.DO, must be created to hold the 


input parameters. DSIN.DO will contain the following variables in the following 


eNS, NPS Si 52, iv 
e The sensor location/parameter(s) matrix, X(NS,NP + 5). 


An example of a input file is shown in Figure 3.2 along with a graph of the 


corresponding sensor coverage areas. Entries in the first line of Figure 3.2 indicate that 


there are two sensors with three parameters each, that the target location distribution 


is BVN(O, 0, 625, 625, 0), and that all sensors have the same parameters. The second 


line is the parameter set for the first sensor. There is no variance in the location of 
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»« SENSOR COVERAGE (SHADED AREAS) 
« TS6T LOS IS DIST BVN(0,0,625,625,0} 
LOCATION DISTRIBUTION (DOTTED CINES) 





Figure 3 ie Example of Input File. DSIN, 
me Graph Of Sensor Coverage. 
either sensor. The first sensor is located at (20,20) and has additional parameters 5, 25, 
and 30. If this input file represents sensors that have a deterministic detection 
functions, then the sensors detect all targets within a distance of 5 and within a circular 
band from 25 to 30 and do not detect targets outside of these bands. Both sensors 
have the same parameters other than location, as indicated by the last element of the 
first line being one, the last line contains only the location of the second sensor. If one 
or more sensors had different parameters (other than location), then the last element in 
the first row would have been zero. Also, the three nonlocation parameters would 
have been specified for all sensors. Although the parameters for each sensor may be 
different, the detection function for each sensor must be the same. For example, the 
program does not allow a mixture of sensors with deterministic and probabilistic 
detection functions. 
2. Interactive Input 
After the program has read the input file, the operator will be prompted to 

specify: 

e Whether a deterministic or a probabilistic detection function will be used. 

e Whether the variance in target location is or 1s not equal to zero. 


e Whether a Monte Carlo simulation or numerical approximation will be used. 


Js 


If the probabilistic detection function is used, the operator may interactively change 
the detection function. If the Monte Carlo simulation is used, the operator will be 


prompted to specify the number of repetitions, NR, in the simulation. 


D. OUTPUT 

The program prints the probability of detection, PD, to the screen of the M100. 
If the computation was a Monte Carlo simulation, the program asks the operator to 
specify an a? for the confidence interval on the probability of detection. Permissible 
a@’s are .10, .05, and .01. The program computes the lower and upper limits of the 
confidence interval based on a normal approximation to the binomial distribution. The 


expression for these limits is shown in Equation 3.3. 
PD + Z gC PD(I-PD)/NRI” (eqn 3.3) 


Ay es ; is the point on a standard normal distribution at which the CDF has a value of 
1l-a/2. 


FE. EXPLANATION OF PROGRAM COMPONENTS 
A complete program listing is at Appendix A. 
l. Initialization/Data Input, Figure 3.3 

Lines 100-115 print the program title and open the input file, DSIN. 

Line 120 reads the first line of DSIN and calculates the variances of the 
marginal distributions of the target location distribution. It also dimensions 
X(NS,5+ NP) and A2(6,6). 

Lines 125-135 read the lines of DSIN that specify parameters for individual 


sensors, S:, 1=1,2,...,.NS. Line 125 reads all parameters for Sy}. If some S; have 


y 
different nonlocation parameters, i.e. F1=0, then line 128 reads in 5+ NP parameters 
for S; where 1=2,....NS. If all sensors have the same nonlocation parameters, i.e. 
F1l=1, then lines 132-135 read the five location parameters for each sensor. They also 
make nonlocation parameters of S;, where 1 = 2,3,4,...,NS, equal to the nonlocation 


parameters of S). 


* 


5a is the probability that the computed confidence level does not contain the true 
probability of detection. 


ae 


CLS: PRINT": PRINT" DETECTION SIMULATION": FORI=1T0400:NEXTI 
‘Input/Initialization 

OPEN'DSIN"FORINPUTAS1 

INPUT#1  NS»NP »S1,S2,RH,F1: V1=S1*S1:V2=S2*S2 

DIMX(NS »54#NP ),A2(656),T1(3 ) 

FORI1=1TO5S+NP: INPUT#1,X(1,I11):NEXTI1 

IFNS=1THEN140 


IFF1=1THEN132 

FORI1=2TONS: FORI2=1TO5+NP : INPUT#1,X(1I1,12):NEXTI2:NEXTI1:GOTO14¢0 
FORI1=2TONS: FORI2=1TO5: INPUT#1 ,X(11,12):NEXTI2: IFNP=OTHEN135 
FORI2=6TOS+NP;: X(11,12)=X(1,12):NEXTI2 

NEXTI1 

RF=SQR(1-RH% 2) 

DFS="EXP(-( (XT-XS 4 241YT-YS )% 2/0 2*X( 2 56% 2)" 





Figure 3.3. Initialization and Data Input Section. 


Line 140 calculates RF, an intermediate value used in later integration 
calculations. | 
Line 150 assigns the right hand side of the Carlton equation to DFS as the 
default equation for calculating the detection function value, DF. 
2. Method Selection Section, Figure 3.4 


‘xX Simulation Selection Section *** 

CLS:PRINT"Is the Detection function:" 

PRINT" 1. Deterministic":PRINT" 2. Probabilistic" 

INPUT“Enter 1 or 2:"3F1 

CLS:PRINT"Are Sensor Locations:" 

PRINT" 1. Always At Aim Point™:PRINT"“ 2. Distributed BVN Around Aim Point" 


INPUT“Enter 1 or 2:"3F2 

IF( F1l=20RF 2=2 JTHENF3=1:GOTO230 

CLS:PRINT':PRINT"Is the Calculation: " 

PRINT" 1=Monte Carlo Simulation":PRINT" 2=Numerical Approximation" 
INPUT“Enter 1 or 2"3F3 

TIMES="00:00:00":IF FL=1THENGOSUB300ELSEGOSUB500 

GOTO200 





Figure 3.4 Method Selection Section. 


Lines 200-250 assign values to: 
e Fl. If Fl=1, then calculations are done for a deterministic sensor. If Fl=2, 


then calculations are done for a probabilistic sensor. 
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e F2. If F2=1, then the sensor is always located at its aiming point, Le. 


67, = a”, =0. If F2= 2, then oa or cae 


e F3. If F3=1, then a Monte Carlo simulation is conducted. If F3=2, then a 


is not equal to 0. 


numerical approximation is calculated. 
Line 230 sets the M100’s clock to zero to time the calculation. It also branches to the 
appropriate subroutine depending upon which type of sensor is specified by F1. 


3. Deterministic Sensors, Lines 300-360 


300 ‘Deterministic Sensor Subroutine 


305 IFF3=1THENGOSUB310ELSEGOSUB350 
306 RETURN 





Figure 3.5 Decision Logic For Deterministic Sensors. 


Based upon the value of F3, lines 300-306 (see Figure 3.5) branch to 
subroutines to conduct the calculations using a Monte Carlo simulation or a numerical 
approximation. 

a. Actual Detection Function, Lines 310-360 

(1) Monte Carlo Simulation, Figure 3.6. Line 315 calls subroutine 900 


which prompts the operator for the number of repetitions, NR. 


"Monte Carlo of Deterministic Sensor 

GCSUB900 

PD=0; FORJ1=1TONR: PRINT. 241,"Repetition: '"3J1:GOSUB600 : FORJZ=1TONS 
IFF2=2THENXS=X(J2,1):YS=X(J2,2)}:GOTO325 

GOSUB612 


T1=SQRE (XS-XT I“ 240 YS-YT )%2) 

IFT1<=X(J2,6 JTHENPD=PD+1:GOT0335 

IFT1>=X( J2,7 JANDT1<=X( J2,8 JTHENPD=PD4+1:GOTO335 
NEXTJ2 

NEXTJ1:PD=PB/NR:GOSUB 950 : RETURN 





Figure 3.6 Monte Carlo Simulation. 


For repetitions one through NR, lines 320-335 generate a target 


location from the BVN distribution. They then check whether that location lies within 
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the circular detection bands around each sensor, S;, l= LQ NS? PD counts: the 
number of repetitions for which the target is in at least one detection band. The 
Monte Carlo simulation functions accurately even if detection bands of various sensors 
overlap because it does not double count if a target is detected by more than one 
sensor. When all repetitions are completed, the counter, PD, is divided by NR to 
produce an estimate of the probability of detection. Finally, output subroutine 950 is 
called. 

(2) Numerical Approximation, Figure 3.7. The volume* under the target 
location distribution is calculated for the detection bands of each S; This volume is 
the sensor’s probability of detection. Overlapping detection bands are not permitted with 
this numerical approximation because the program would double count the overlapping 
volumes. Subroutine 1200 conducts the integration. Probability of detection is the 


volume under the target location distribution that is within the coverage area of any S.. 


"Numeric/Deterministic Subroutine 
PD=0: FORJZ2=1TONS: H=X(J2,6):GOSUB1200:PD=PD+IN 


H=X(J2,8):GOSUB1200: PD=PD+IN 
H=X(J2,7):GOSUB1200:PD=PD-IN:NEXTJ2 
GOSUB 950: RETURN 





Figure 3.7 Numerical Approximation. 


4. The Probabilistic Sensor, Figure 3.8 

Lines 502-503 print a header and call subroutine 1300 which permits the 
detection function to be modified. Line 503 also calls the subroutine in which the 
operator specifies the number of Monte Carlo repetitions. 

For each repetition, lines 520-530 generate a target location from the BVN 
target location distribution. Then, for each S; they calculate the value of the sensor's 
detection function using subroutine 1410 and generate a uniform random number 
between zero and one. If that random number is less than or equal to S;’s detection 


function value at that location, then S; detected the target and PD is augmented by 


4This volume, although it is computed and described as a volume in this section, 
is not a true volume. This is because although the X and Y_axes of the BVN are 
distances, the Z axis is a probability, i.e. dimensionless. Therefore the result is really 
an area, not a volume. 


zo 


‘Probabilistic Detection Function 

CLS:PRINT' Default Detection Function Is Carleton." 
GOSUB1300 : GOSUB900 

PD=0: FORJ1=1TONR: PRINT .241,"Repetition: "3J1:GOSUB600: FORJZ2=1TONS 
IFF2=2THENXS=X(J2,1): YS=X(J2,2):GOTO523 


GOSUB612 

GOSUB1410: IFRNS(1 )<=DFTHENPD=PD +1 :GOTO526 
NEXTJ2 

NEXTJ1:PD=PD/NR 

GOSUB 950: RETURN 





Figure 3.8 The Probabilistic Detection Function. 


one. If one sensor detects the target, the program moves directly to the next 
repetition. Therefore, the Monte Carlo simulation functions accurately even if detection 
bands of various sensors overlap because it does not double count a target that is detected 
by more than one sensor. When all repetitions are completed, the counter, PD, 1s 
divided by NR to produce an estimate of the probability of detection. Finally, output 
subroutine 950 is called. 

5. Generating A BVN Random Variable, Lines 600-606 


*xx¥*¥Generate BYVN RV*x*x* 
U1=RND(1 ):U2=RND( 1): TE=SQR( -~-2*LOG( U1 )) 
XT=TEXCOS(6. 2831853 *U2 ): YT=RH¥XT +RFXTEXSIN( 6. 2831853*UZ ) 


XT=XT¥S1:YT=YT#S2: RETURN , 

U1=RND(1):U2=RND( 1): TE=SQR( -2*LOG( U1) ) 

XS=TEXCOS( 6.2831853U2 ): YS=X( U2 »5 )*XT+01-X0 52 5 *2)% . SRTEXSIN( 6. 2831853%U2 ) 
XS=XCI2 91 4+XS¥*X0 253 )2 YS=X( S222 V4YS*X( JZ >): RETURN 





Figure 3.9 Generating A BVN Random Variable. 


Lines 602-606 generate the target location. Lines 602-604 compute both 
components of a BVN(0, 0, 1, 1, 0) random variable using two independent uniform 
(0,1) random numbers generated by the M100’s RND(1) function. Equations 3.4, and 
3.5, as described in [Ref. 4:page 953], form the basis for lines 602 and 604. Equations 


3.4 generate two independent normal(0,1) random variables, X,’ and X4’. 


X)/ = (-2ln U,)> cos(2nU) (eqn 3.4) 
Xo’ = (-2In Uy)? sin(2nU,) 
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Equations 3.5 convert these two independent normal random variables into the 


components of a BVN(0, 0, 1, 1, p) distribution, Xx] and X>. 


ag, ae Legit) 3.9) 
X= Oy) ORS 


Line 606 scales the components of the BVN(0, 0, 1, 1, Px y) by S1 and S2 to 
produce the components of the BVN(O, 0, ae ace, Py y) target location distribution 
and then returns to the calling program. 

Lines 612-616 generate a sensor location using the same algorithm that was 
used to generate the target location. However, in addition to being scaled by 6, and 
G,,, sensor locations must also be displaced by t,, and HL, 


6. Entering The Number Of Repetitions, Figure 3.10 


900 C!S: INPUT"Enter number of repetitions for Monte Carlo Simulation:"3NR 


905 RETURN 
910 INPUT"Hit ENTER to Continue"3Z1:RETURN 





Figure 3.10 Entering The Number Of Repetitions. 


Line 900 prompts the operator to interactively specify the number of 
repetitions, NR, for Monte Carlo simulations. 

Line 910 is a subroutine which stops the program while the operator views 
output to the screen. 

7. Output Section, Figure 3.11 

All printing is to the screen of the M100. Lines 951-952 play a tune to notify 
the operator that the calculation is finished. Line 953 also prints the time required to 
do the calculation and branches around the selection of a for the confidence interval if 
a numerical approximation was used. Lines 954-959 prompt the operator to select .1, 
.05, or .01 as @=AL, the probability that the true probability of detection is not in the 
confidence interval. Once AL is selected, it is reassigned the value of the standard 
normal at Z 1.4/2: Line 960 prints the point estimate of Py. Line 962 prints a short 


reminder that there is no confidence interval with numerical approximations. Line 966 


Ze 


‘Print output 

SOUN01567 , 10: SOUN01244,10:SOUN01046 »10:SOUNO783 ,20 

SOUN01046 ,10:SOUNO0783 540 

CLS:PRINT'':PRINT"Calculation Time (HH/MM/SS) = "$3TIME$:IFF3=2THEN960 
PRINT"Select Alpha for Confidence interval:" 

INPUT" Choices = .1, .05, .01:"3AL 

IFAL=.1THENAL=1.645:GOT0O960 

IFAL=.O5THENAL=1.96:GOTO960 

IFAL=.O1THENAL=2.575:GOTO960 

GOTO954 

PRINT"*** Estimate of P(Oetection) = '"3:PRINTUSING" #.#####"'53PO0O 
IFF3=1THEN965 

PRINT"'No Confidence Interval For Numerical Approximations” 

GOTO0970 

PRINT"Confidence Interval: ("5 
TE=AL*SQR(PD*(1-P0 J/NR): LL=PO-TE :UL=PO4TE : IFUL>1LTHENUL=1 
IFLL<OTHENLL=0 : 
PRINTUSING"## . #####"5LL3UL3:PRINT" )":GOSUB910 

"Confetti Approximation 

PRINT": INPUT"Confetti approximation? O=No, 1=Yes:"3Z29: IFZ9=OTHENRETURN 
CLS: INPUT"Enter TOTAL lethal area for ALL sensors in the pattern:"3NA 
TE=NA/(6.283185*S1*S2 ): TE=1-(1+SQR( 2*TE ) JXEXP( -SQR( 2*TE )) 
PRINT"*xConfetti Approximation = ";TE:GOSUB910:RETURN 





Figure 3.11) Output Secon, 


calculates the confidence interval according to Equation 3.3. Lines 966-967 ensure that 
the confidence limits are between zero and one. Lines 965 and 968 print the confidence 
interval. 

Lines 970-977 approximate P, with the “confetti approximation” described in 
[Ref. S:pp. 14-16]. This approximates P, by distributing the total lethal area of the 
entire group of sensors over an ellipse as described in Reference 5. This approximation 
is calculated by Equation 3.6 where ¢ = (Total Lethal Area)/(206,6y). 
P ge ee (20)>) (26) (eqn 3.6) 


Line 972 prompts the operator to indicate whether a confetti approximation is 
desired and ends the output routine if it is not. Line 974 prompts the operator to enter 
the total lethal area for all sensors combined. Line 976 computes § and then P,. Line 
977 prints P, and ends the output subroutine. 

8. Numerical Integration, Figure 3.12 

This subroutine is a tailored version of the Romberg integration subroutine 

documented in Chapter 8. It integrates the target location distribution subject to 


circular limits of integration. 
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‘Numerical Integration Subroutine 

D1=6 .2831853*S1*S2*RF 

TL=.001 

CLS: PRINT: PRINT" t!Calculating An Integral! !":PRINT"" 
N=2:GOSUB1293 :DY=(YU-YL)/2 

FORJ9=1T06 : DY=DY/2 : N=N*2 

Y=YU: GOSUB12 96 : GOSUB1280 : A2(J9,1)=TS*DX 

Y=YL:GOSUB12 96 : GOSUB1280: A2(J9,>1 J=A2(J9,1}+TS*¥DX 
FORJ8=2TON: Y=Y+DY : GOSUB12 96 : GOSUB1280 

A2(J9,1 3J=A2(J9,1)4#2*TS¥DX: NEXTUS 

A2(J9,1)=A2(J9,1 )¥DY/2 

IFJ9=1THENNEXTU9 

FORJ8=1TOJ9-1 

A2(J9,J8+1 J=A2(J9,J8 4 (A209, J8 J-A2( S9-1 8 ) 3/14 J8-1 } } = NEXTUS 
T1=A2(0J59,J9 )-A2( 59 J9-1 Js IFSGN(T1 IXTL-TL>OTHENNEXTJ9ELSE1266 
PRINT"Tolerance of"3TL3"not met after five extrapolations" 
IN=A2(J9,J9}: RETURN . 

FORJ7=1T06: FORJ6=1TOJ7: PRINTUSING''S# . #34" 5A2( U7 ,J6 33 

NEXTJ6: PRINT": NEXTJ7: INPUTZ9: RETURN 

REM Trapezoidal Rule Sum 

X=XU: GOSUB1286 : TS=F :X=XL: GOSUB1286: TS=TS+F 

FORJ5=2TON-1: X=XK+DX; GOSUB1286: TS=TS+F :NEXTJ5: RETURN 

‘F(X»,Y} to be integrated: 

F=X* 2/V1-2*RH*¥X*Y/S1/S2+Y% 2/V2 

F=(EXP(-F/2/RF* 2} }/D1:RETURN 

‘Limits of Integration: 

YU=X(J2,2 )4H: YL=X(J2>,2 }-H: RETURN 

T3=SQR(H* 2-(Y-X( 5252 ))% 2): XU=X( S251 473: XL=EX( J251 9-73: DX=(XU-XL I/N 
RETURN 





Figure 3.12 Numerical Integration Subroutine. 


9. Changing The Detection Function, Figure 3.13 


1300 PRINT" -Detection Fn (DF} in terms of XT, YT," 
1302 PRINT" and Parameters XS, YS, and X(J2,6})>.... ,X(J2,5+NP):" 


1304 PRINT’ ** DF = ";DF$ 
1306 PRINT"Hit ENTER For No Change or Enter New...":INPUT" DF = "3DF$ 
1307 RETURN 





Figure 3.13 Changing The Detection Function. 


This section contains a subroutine which permits the operator to interactively 
change the detection function, lines 1300-1306., These changes would be made if the 
Operator wanted to use a probabilistic detection function other than a Carlton 


function. In both equations XT and YT represent the two components of the target 
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location. XS and YS represent the location of Sy5 in the same coordinate system as 
XT and YT. An example of a probabilistic detection function that is not Carlton 


might be an exponential function as specified in Equation 3.7. 
Pa he! where r = distance from sensor to target. (eqn 3.7) 


The entry to be made when prompted by line 1306 would be 
X(J2,6)"*EXP(-X(J2,6)*SQR(U(XT-XS)*2 + (YT-YS)*2)) 


where X for sensor J2 is X(J2,6). In general, X(J2,6),...,.4(J2,5-+ NP) represent NP 
other parameters of the selected detection function in addition to the five parameters of 
the BVN sensor location distribution. The formula for the current detection function 1s 
displayed by line 1304. The operator may either change it as desired or hit ENTER to 
leave the current formula unchanged. 


10. Formula Tokenization Section, Figure 3.14 


1400 'ToKenize OF 
1410 BS="DF="4+DFS+CHRS(0O } 
1450 'ToKenize/execute BS 


1451 BO=VARPTR(BS ):B1=PEEK(BO+] }+256*XPEEK(BO+2 }:CALL1606,0,B1 
1455 CALL2499,0,63105: RETURN 





Figure 3.14 Formula Tokenization Section. 


The right hand side of the detection function equation is stored as a string 
variable, DF$. This section converts DFS into an executeable BASIC assignment 
statement and executes that statement. A detailed explanation of this subroutine is 


found in the Formula Tokenization Section in Chapter 2. 


F. SAMPLE PROBLEMS 
Examples 1-6 which follow have the following characteristics in common. 
e They use a target location distribution that is BVN(O, 0, 252, 252, QO), except for 
Example 6 in which Priya A: 
e All measurements are in kilometers and are in the coordinate system of the target 
location distribution. 


e All confidence intervals are calculated with @=.05. 
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¢ All Monte Carlo simulations were done with 5,000 repetitions. 
e <All sensors, deterministic or probabilistic, have a lethal radius of twenty and 
therefore a lethal area of 400n. 
1. Example | 
This example is of a single cookie-cutter sensor located at (0,0) with a lethal 
radius of 20. The input file and diagram of sensor coverage is at Figure 3.15. The 


closed form solution may be calculated using Equation 3.8. 


Dep 2192952 
Py = l-e7 120° = 1. 20°(2"25°) = 97385 (eqn 3.8) 


m SENSOR COVERAGE (SHADED AREAS) 
mn TGT LOC IS OIST BVN{0,0,525,625,0 
LOCATION DISTRIBUTION (DOTTED LINES) 





Figure 3.15 Input File And Sensor Coverage Diagram For Example 1. 


The probabilities of detection computed by the various computational 
techniques and calculation tmes were as follows. 
Seelosed form: .2/385 
e Numerical integration: .27252; 3 minutes, 24 seconds 
e Monte Carlo simulation: .27640 + .0124; 4 hours, 22 minutes. 
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2. Example 2 
This example is of a single Carlton sensor located at (0,0) with a lethal area 
equivalent to the cookie-cutter sensor in Example 1; 1e., b= 20/(2°>) = 14.14214. The 
input file and diagram of sensor coverage is shown in [Ref. 5:page 5]. The closed form 
solution is at Equation 3.9. Because the sensor is located at the origin and o, = o 


X y; 
Equation 3.9 simplifies to Equation 3.10 for this example. 


1125 2500 Carlton Sensor (b=14.1421) Located At (0,0) 
00000 14.14214¢ 
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Figure 3.16 Input File And Sensor Coverage Diagram For Example 2. 


The probabilities of detection computed by the various computational 
techniques and calculation times were as follows. 
¢ Closed form: .242424 
e Numerical Integration: None 


© Monte Carlo simulation: .2452 + .0119; 3 hours, 44 minutes. 


Py og f e7 5(84, + Oy) 
where Y = ba/(b- +66 coe 


5, = nw /(b? + 07,), and 6,=pn2,/(b2+02,). 


(eqn 3.9) 


SZ 


Pa = ba(b? + 6%) = 14.1422/(14.1422 + 252) = .242424 (eqn 3.10) 


3. Example 3 
This example is for a single cookie-cutter sensor offset at (/0,/0) with a lethal 
radius of 20. The input file and diagram of sensor coverage is shown in Figure 3.17. 


There 1s no closed form solution for offset cookie-cutter sensors. 


13 25 2500 


10 100002000 » SENSOR COVERASE (SHADED AREAS) 
x TGT LOC 15 OIST BYN(0,0,525,625,0) 
LOCATION DISTRIBUTION (COTTED LINES) 
[ = | 


2 





Figure 3.17 Input File And Sensor Coverage Diagram For Example 3. 


The probabilities of detection computed by the various coimputational 
techniques and calculation times were as follows. 
¢ Closed form: None. 
e Numerical integration: .2400; 2 minutes, 51 seconds. 
e Monte Carlo simulation: .24160 + .01187; 2 hours, 57 minutes. 
4. Example 4 
This example is for a single Carlton sensor offset at (10,10) with a lethal area 
equivalent to the cookie-cutter sensor in Example 3, 1.e. b=20/(2>)= 14.14214 The 
input file and diagram of sensor coverage is shown in Figure 3.18. The equation for 


the closed form solution is shown in Equation 3.9. 
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Figure 3.18 Input File And Sensor Coverage Diagram For Example 4. 


The probabilities of detection computed by the various computational 
techniques and calculation times were as follows. 
¢ Closed form: .21475 
e Numerical integration: None 
e Monte Carlo simulation: .2138 + .0114; 3 hours, 15 minutes. 
5. Example 5 
This example is for a single convergence zone sensor located (0,0) with a lethal 
area equivalent to the sensors in Examples 1-4. The sensor detects all targets at ranges 
less than four, detects no targets between four and 30, detects all targets between 30 
and 35.83, and detects no targets beyond 35.83. The input file and sensor coverage 
diagram is shown in Figure 3.19. The closed form solution is found by using Equation 
3.8 to calculate Py in circles of radu 4, 30, and 35.83. Then Pg = Pg(4) - Pg(30) + 
P 3(35.83). 
The probabilities of detection computed by the various computational 
techniques and calculation times were as follows. 
¢ Closed form: .141463 
e Numerical integration: .14128 
e Monte Carlo simulation: .1416 + .0097; 2 hours, 59 minutes. 
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0 35.83295 » SENSOR COVERAGE (SHADED PREPS) 
« TST LOC 1S DIST BVN(9,0,625,6 
LOCATION DISTRIBUTION BOrteD: (INES) 





Figure 3.19 Input File And Sensor Coverage Diagram For Example 5. 


6. Example 6 
This example is for four convergence zone sensors with mean locations at 
(10,10), (-10,-10), (10.-10), and (-10,10). These sensors have convergence zones equal 
to that of the the sensor in Example 5. Sensor locations are distributed BVN with pL, 
and H, as indicated above and 6,,=6,=3, and Puye fee ie inp Miler iSear) rE ieune 


3.20. There is no closed form solution for this example. 


So 25 25 75> 2 
10 10 33 .7 4 30 35.83295 





Ricine J c0sernputiaile for Example 6. 


The probabilities of detection computed by the various computational 
techniques and calculation times were as follows. 
e Closed form: None 


¢ Numerical integration: None 
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e Monte Carlo simulation: .4570 + .01381; 2 hours, 53 minutes. 
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IV. KALMAN FILTER PROGRAM 


A. GENERAL 

This program successively updates an estimate of the state of a system, [t, based 
upon a series of measurements, ¢. Vectors ft and € need not have the same dimensions. 
The program makes provision for changes in the system state between measurements in 
accordance with a linear system model. Covariances between elements of HM and 
between elements of ¢ are also accounted for by the program. 

The purpose of this chapter is to describe an implementation of a Kalman filter 
on the M100. A fuller explanation of the mathematical background behind Kalman 
filters may be found in References 6 and 7. The notation in the following program is 
similar to the notation used in Reference 6 to facilitate comparison of the 
mathematical theory and computer implementation. 

The state of the system, X, is assumed to be a multivariate normal random 
variable, X ~ N(p,2), with system noise, W ~ (jl,,,Q), and measurement noise, V ~ 
(H),R). @ is the matrix which models the linear change in X between measurements. 
H is the matrix which shows the linear relationship between the measurement and the 
system State, 1.e. how the measurement depends on the system state. 

The Kalman filter recursively updates t by repeating two steps, measurement and 
movement. The measurement step calculates the Kalman gain, enters a new 
measurement, and updates #t and X based on that measurement. The movement step 


updates ft and LX based upon the system model. 


B. EXPLANATION OF VARIABLES 
e Bl is the intermediate matrix for inversion of C2. 


¢ Cl(MD,MD) and C2(MD,MD) are the matrices which hold intermediate results 
of matrix calculations. 


e H is the matrix showing the relationship between measurements and the system 
State: 


oan i |? dneelo are loop counters: 

¢ K(NX,NZ) is the Kalman gain matrix. 

e MD is the maximum of NX and NZ. 

¢ MU(NX) is the current estimate of the system state. 
¢ MV(NZ) is the mean of measurement noise. 


¢ MW(NX) is the mean of system noise. 


SY 


input. 


NX is the number of elements in the system state vector. 
NZ is the number of elements in the measurement vector. 


PH(NX,NX) is the matrix, @, modeling the linear change in the system state 
between measurements. 


Q(NX,NX) is the covariance matrix of system noise. 
R(NZ,NZ) is the covariance matrix of measurement noise. 
SG(NX,NX) is the covariance matrix of the system state. 


Z(NZ) is the vector holding the current measurement, ¢. 


INPUT 

The operator must create a RAM file, KALIN.DO, which contains the program 
KALIN.DO must contain the following parameters and matrices in order. 

NX and NZ, the number of elements in wt and € respectively. 

Matrix PH(NX,NX), the system model matrix, @. 

Vector MW(NX), the mean of system noise, [M,,. 

Matrix Q(NX,NX), the covariance matrix of system noise. 

Matrix H(NZ,NX), the matrix showing the relationship between pi and C. 

Vector MV(NZ), the mean of measurement noise, }1,. 

Matrix R(NZ,NZ), the covariance matrix of measurement noise. 

Vector MU(NX), the initial estimate of the system state. 


Matrix SENS, NX), the initial estimate of the covariance matrix of the system 
state, 


After the initial estimate of the system state is entered from the input file, 


KALIN.DO, the program will prompt the operator for measurements and changes to 


the H matrix to be entered from the keyboard. 


D. 


OUTPUT 
All output goes to the screen of the M100. After the measurement step the 


program displays the updated Kalman gain matrix, the estimated system state, h+, 


and the covariance matrix, 2+. After the movement step the program displays the 


estimate of the system state, [fi-, and covariance matrix, L-, just prior to the next 


Measuremecnl 


| 


EXPLANATION OF PROGRAM 
A complete program listing is located at Appendix B. 
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1. Initialization/Input, Figure 4.1 


CLS: PRINT" KALMAN FILTER": PRINT" Input Data Being Read" 
OPEN" KALIN"FORINPUTAS1 : ONERRORGOTO9900 

INPUT#1 NX »NZ: IFNX<NZTHENMD=NZELSEMD=NX 

DIMPH( NX »NX ) »MWIONX ) >QONX NX ) >HONZ »NX ) »DMVONZ ) »RONZ NZ) »>MUCNX ) >SGONX NX ) 
DIMC1(MD MD) ,C2(MD MD ) KC NX>NZ ) ,B1(NZ4+1)NZ2 ) 

FORIL=1TONX: FORI2=1TONX: INPUT#1,PH(I1,I2):NEXTI2:NEXTI1 
FORIL=1TONX: INPUT#1 »MW(I1):NEXTI1 


FORT1=1TONX: FORI2=1TONX: INPUT#1,Q(I1,12):NEXTI2Z:NEXTI1 

FORIL=1TONZ: FORT 2=1TONX: INPUT#1,H(I1,12):NEXTI2:NEXTI1 

FORI1=1TONZ: INPUT#1 »MV(I1):NEXTI1 

FORIL=1TONZ: FORI2=1TONZ: INPUT#1 >R(I1,I2):NEXTI2:NEXTI1 
FORIL=1TONX: INPUT#1 ,MU(I1):NEXTI1 

FORIL=1TONX: FORIZ2=1TONX: INPUT#1 ,»SG(I1,1I2):NEXTI2:NEXTI1 
CC=0:CLS:PRINT"Initial SG As Input Check: ":GOSUB532 





Figure 4.1 Initialization and Input Section. 


Line 110 opens the input file, KALIN.DO and branches the program to line 
9900 if an error occures. Line 120 enters NX and NZ and calculates MD. Lines 125 
and 126 dimension the matrices in the program. Lines 130-144 enter the matrices from 
KALIN.DO as described in the input section above. Line 145 initializes the 
measurement counter, CC, and prints the last input matrix, SG, as a check on input 
accuracy. 

2. Measurement Block, Lines 150-387 

Measurement block matrix equations for updating K, pm, and 2 are listed in 

Equations 4.1, 4.2, 4.3. 


K = SHYHEH' + Ry! (eqn 4.1) 
ee KZ = ht, - Hp) (eqn 4.2) 
pe (lo (eqn 4.3) 
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G. Entering A New H Matrix, Figure 4.2 


CLS:PRINT" HHHXHXHXMEASUREMENT BLOCK 2¢2¢3¢2¢2%"" 
PRINT"Current 4H '":GOSUB540 
INPUT"Enter New H ? 1=Yes, O=No:"329: IFZ9=OTHEN170 


‘Enter A New H 

FORI1=1TONZ: FORI2=1TONX 

PRINT"Enter Row"3I1;", Column"31I23;"0Of H :%3 
INPUTH(I1,I2):NEXTI2:PRINT":NEXTI1 





Figure 4.2 Entering a New H Matrix. 


Some Kalman filter problems require a different H matrix for each 
measurement. This section permits such a matrix to be entered. Line 160 prints a 
header and calls the subroutine which prints the current H matrix. Line 162 asks the 
operator whether a new H matrix is required and branches the program appropriately. 
Lines 167-169 prompt the operator to fill the new H matrix row by row. 


b. Kalman Gain Calculation, Figure 4.3 


"CALC KALMAN GAIN 

"MULT SG H t», INTO Cl 

FORI1=1TONX: FORI2=1TONZ:C1(1I1,12 )=0: FORIZ=1TONX 

C1(11,12 )=(SG(I1,13 )*H(I2,13 ))+C1(I1,I2):NEXTI3: NEXTI2:NEXTI1 
"MULT H SG Ht » INTO C2 

FORI1=1TONZ: FORI2=1TONZ:C2(I1,I2 )=0: FORI3=1TONX 

C2(1I1,12 J=(H(I1,13)*C1(1I3,12))+C2(1I1,12):NEXTI3:NEXTI2:NEXTI1 
"ADD R_ INTO C2 

FORI1=1TONZ: FORI2=1TONZ:C2(1I1,12 )=C2(1I1,12 )+R(1I1,12) 
NEXTI2:NEXTI1 

"INVERT C2 

GOSUB 9800 

"MULT Cl C2 INTO K 

FORI1=1TONX: FORI2=1TONZ:K(1I1,I2 )=0: FORI3=1TONZ 
K(I1,I2)=(C1(I1,13 )*C2(I3,12) +K(I1,12):NEXTI3:NEXTI2:NEXTI1 





Figure 4.3. Kalman Gain Calculation. 


This section updates the Kalman gain matrix, K, in accordance with 
Equation 4.1. Lines 171-174 multiply £ by H' and place the result in matrix Cl. 
Lines 180-184 multiply H by ZH! and place the result in matrix C2. Lines 200-203 add 
R to (HEH) and place the result in C2. Line 215 calls the subroutine which inverts 
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(HEH + R). Lines 220-224 multiply ZH" by (HEH' + Ry, producing an updated 
Kalman gain matrix, K. 


c. Enter A Measurement and Update yp, Figure 4.4 


"*¥*X*X*XXUPDATE MU- TO MU +2 

"MULT H MU- INTO Cl 
FORI1=1TONZ:C1(I1,1 )=0: FORIZ=1TONX 

CX.(I1,1 JS0CH(I1,I3 )XMU(I3 ) 4+01(I1,1):NEXTI3:NEXTI1 
‘ADD MV + H MU- 
FORII=1TONZ:C1(I1,1)=C1(I1,1)4+MV(I1):NEXTI1 

"INPUT A NEW MEASUREMENT 
CC=CC+1:CLS:PRINT"Measurement #'"3CC3";" 
FORI1=1TONZ:PRINT"Enter Element"3I13;"0f Measurement: ''3 
INPUTZ(I1):NEXTI1 

‘SUBTRACT C1 FROM Z, INTO Cl 

FORII=1TONZ:C1(I1,1 )=Z(I1)-C1(I1,1):NEXTI1 

"MULT K C1 INTO C2 

FORIL=1TONX:C2(1I1,1)=0: FORIZ=1TONZ 

C2(I1,1 J=(K(I1,I3 )*¥C1(I3,1))+02(I1,1):NEXTI3Z:NEXTI1 
"ADD C2 + MU- TO UPDATE TO MU+ 
FORI1=1TONX:MU(I1)=C2(I1,1)+#MU(I1):NEXTI1 





Figure 4.4 Enter Measurement And Update The Estimate of ph. 


This section enters a new measurement, Z, from the keyboard and updates 
the estimate of pt in accordance with Equation 4.2. Lines 251-254 multiply H by pt and 
place the result in Cl. Line 262 adds ft, + Hp and places the result in Cl. Lines 
272-274 increment the measurement counter and allow the operator to enter the new 
measurement, Z. Line 282 subtracts (ft, + Ht) from Z, and places the result in Cl. 
Lines 292-294 multiply the Kalman gain, K, by (Z - ph, - Hp), and place the result in 
C2. Line 302 adds pt to K(Z - Wt, - Hpt), producing the revised estimate of }t. 

d. Updating &, Figure 4.5 

This section updates the estimate of 2 using Equation 4.3. Lines 322-326 
multiply the Kalman gain, K, by H, subtract the result from the identity matrix, I, and 
place the result in Cl. Lines 352-354 multiply (I - KH) by 2 and place the result into 
C2. This result, the updated 2, is then copied into SG by line 362. 

e. Printing The Updated K, jt, And &, Figure 4.6 

Lines 375-377 print a header and call the printing subroutine for the 

updated Kalman gain. Lines 380-382 print a header and call the printing subroutine 


for the updated estimate of the system state, ft. Lines 385-387 print a header and call 
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‘MULT K H & SUBTR FROM I ,», PUT IN Cl 

FORI1=1TONX: FORI2=1TONX:C1(I1,I2 )J=0: FORIZS=1TONZ 

C1(I1,122 J=(K(I1,13 J¥XH(I3,12))4+C1(11,102 ):NEXTI3:C1l(I1,I2)=-C1(I1,1I2) 
NEXTI2:NEXTI1 

FORII=1TONX:C1(1I1,1T1)=1+C1(I1,I1):NEXTI1 


‘MULT LAST RESULT BY SG , INTO C2 

FORI1=1TONX: FORI2=1TONX:C2(1I1,I12 J=0: FORIZ=1TONX 
C2(I1,12)=(C1(I1,13 J*¥SG(1I3,12) +C20I1,12 )sNEXTIS:NEXTI2:NEXTI1 
"PUT C2 INTO SG 

FORI1=1TONX: FORT2=1TONX:SG(I1,I2 )=C2(I1,I2):NEXTI2:NEXTI1 





Figure 4.5 Updating The System Covariance Matrix, 2. 


CLS:PRINT"Kalman Gain, K(i,j) After" 
PRINT''Measurement #"3;CC:GOSUB510 
CLS:PRINT"Estimate Of System State, MU(i)+ After" 


PRINT"Measurement #''5;CC:GOSUB520 
CLS:PRINT"Estimate Of Covar, SG(1,j3)+ After" 
PRINT'Measurement #"3CC:GOSUB530 





Figure 4.6 Printing The Updated Kalman Gain, pi, And &. 


the printing subroutine for the updated estimate of covariance matrix of the system 


state, 2. 
3. Movement Block, Lines 400-490 


Movement block matrix equations for updating wf and 2X are listed in 


Equations 4.4, and 4.5. 


w= out wy (eqn 4.4) 


x= oxo’ + Q (eqn 4.5) 
a. Updating The Estimate of jt, Figure 4.7 


Lines 422-424 multiply @ by wt and place the result in Cl. Line 432 adds 


Ht... tO PH, producing the updated estimate of ft just before measurement CC + 1. 
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CLS: PRINT "3x33 MOVEMENT mea 
"Update MU(CC)+ to MU(CC+1)- 
"MULT PH MU , PUT IN Cl 


FORI1=1TONX:C1(1I1,1)=0: FORIZ=1TONX 
C1(I1,1)=(PH(I1,I3 )XMUCI3))+C1(I31,1):NEXTI3Z:NEXTI1 
"ADD Cl+ MN , INTO MU 
FORIL=1TONX: MUCI1 J=C10I1,1)4+MNCI1):NEXTI1 





Figure 4.7 Updating The Estimate of HL. 


b. Updating 2, Figure 4.8 


"*¥XUPDATE SG *x 

"MULT PH SG » INTO Cl 

FORI1=1TONX: FORI2=1TONX:C1(I1,I2)=0: FORIZ=1TONX 

C1(I1,I2 )=(PH(I1,I3 JXSG(I3Z,12)34+C1(I1,12):NEXTIZ:NEXTI2:NEXTI1 
"MULT Cl PH t, INTO C2 

FORII=1TONX: FORI2=1TONX:C2(I1,1I2 )=0: FORIZ=1TONX 

C2(IT1,I2 )=(C1(I1,13 JXPH(I2,13))+C20I1,12) sNEXTIZ:NEXTIZ2:NEXTI1 
"ADD C2 + Q = SG 

FORIL=1TONX: FORIZ=1TONX:SG(I1,12 )=C2(I1,12)+Q(I1,I2) 
NEXTI2:NEXTIL 





Figure 4.8 Updating & 


Lines 452-454 multiply @ by 2 and place the result in Cl. Lines 462-464 
multiply @= by o! and place the result in C2. Lines 472-473 add Q to oZq!, 
producing the updated estimate for X& just before measurement CC + 1. 

c. Printing The Updated p And X, Figure 4.9 


480 PRINT"Estimate Of System State, MU(I)-" 
482 PRINT'Before Measurement #';CC+1:GOSUB520 


485 CLS:PRINT"Estimate Of Covar, SG(I,J)J- Before" 
487 PRINT'’Measurement #''3CC+1:GOSUB530 
490 GOTO160 





Figure 4.9 Printing The Updated p And &. 
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Lines 480-482 print a header and call the printing subroutine for the 
updated estimate of the system state, pt. Lines 485-487 print a header and call the 
printing subroutine for the updated estimate of the system state covariance matrix, L. 
Line 490 branches the program back to the beginning of the measurement block. 

4. Printing Subroutines, Figure 4.10 
Lines 500-554 contain subroutines which print the Kalman gain matrix, pt, 2, 


the H matmx, or the C2 matm: 


‘PRINTING SUBROUTINES 

‘PRINT KALMAN GAIN, K : 

FORI1=1TONX: FORIZ2=1TONZ: PRINTUSING"##H#H##H. FH" 5SK(0I1,I2)3:NEXTI2 
PRINT" :NEXTI1L: INPUT"Hit ENTER To Continue: "3Z29:RETURN 

‘PRINT MU 

FORI1=1TONX: PRINTUSING'###H###. HH" 5MUCI1)3:NEXTI1L: PRINT" 
INPUT"Hit ENTER To Continue:"5;29:RETURN 

‘PRINT COVAR MATRIX, SG 

FORI1=1TONX: FORI2=1TONX: PRINTUSING'#####. ##"5SG(I1,1I2)3:NEXTI2 
PRINT’ :NEXTI1: INPUT'Hit ENTER To Continue:"3Z9:RETURN 

‘PRINT H 

FORI1=1TONZ: FORI2=1TONX: PRINTUSING''H####H#H. HH" 53HCI1L,I2)3:NEXTI2 
PRINT": NEXTI1:RETURN 

PRINT" C2 MATRIX:" 

FORI1=1TOA: FORI2=1TOB: PRINTUSING######.##"5;C2(1I1,12)3:NEXTI2 
PRINT’ sNEXTI1L: INPUT"Hit ENTER To Continue: '"'3Z9:RETURN 





Figure 4.10 Printing Subroutines. 


5. Inversion Subroutine For C2, Figure 4.11 
This subroutine is essentially the same as the matrix inversion subroutine in 
the matrix algebra program discussed in Appendix E. It has been abbreviated to invert 
only matrix C2 instead of inverting several matrices of various dimensions as in 
Appendix E. This subroutine also does not calculate the determinant of C2 to test for 
invertability. If C2 is not invertable, an division by zero error will occur and the 
program will branch to the error identification section. A detailed explanation of how 
the matrix inversion subroutine functions is located in Appendix E. 
6. Error Identification, Figure 4.12 
Line 9900 prints a message indicating that C2 is not invertable. Line 9900 1s 
based upon the assumption that a division by zero error in the inversion subroutine 
means that C2 is not invertable. Line 9905 prints the error code and line number of 


other errors. See page 217 of Reference | for an explanation of error codes. 
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"INVERT C2 

FORT1=1TONZ: FORI2=1TONZ:B1( 11,12 )=C2(I1,I12):NEXTI2:NEXTI1 
FORIL=NZ+1T0O2*NZ: FORI 2=1TONZ 

IFIL=I2+NZTHENB1(I2,I1 )=1LELSEB1(I2,I1 )=0 

NEATI2:NEXTI1 

FORI1=1TONZ 
ML=1/B1(1I1,I1):FORI3=1TO2*NZ:B1(I1,I3)=B1(I1,I3 )*ML:NEXTI3 
IFIL=NZTHEN9865 

FORI2=I1+1TONZ: IFB1(12,I1 )=OTHEN9860 

ML=-B1(1I2,I1) 
FORIS=I1TO2*NZ:B1(1I2,13 )=Bl(12,I3 )+#(ML*B1(I1,I3)):NEXTI3Z 
NEXTI2:NEXTI1 

FORI1L=NZTOZSTEP-1 

FORI2=I1-1TO1LSTEP-1:IFB1(1I2,I1 )=OTHEN9885 

ML=-B1(I2,I1) 
FORI3=17T02*NZ:B1(1I2,I3)=Bl(12,13 )+#(ML*B1(1I1,I3 )):NEXTI3 
NEATI2:NEATI1 

FORI1=1TONZ: FORI2=1TONZ 
C2(I2,I1)=Bl(I2,I1+NZ):NEXTI2:NEXTI1 

MI=1:RETURN 





Figure 4.11 Inversion Subroutine For C2. 


9900 IFERL>9800ANDERR=11THENPRINT"!!!'ERROR: C2 Is Not Invertable!!!":END 


9905 PRINT"Error Code"3ERR3"In Line”"3ERL:END 





Figure 4.12 Error Identification. 
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V. MODELS OF COMBAT USING LANCHESTER EQUATIONS 


A. GENERAL 

This program is an example of a time step force attrition simulation using 
Lanchester equations. The scenario is that there are two sides, refered to hereafter as 
the attacker and the defender.> The battle may be broken into phases if some model 
parameters change during the course of the battle. Each side has a fixed number of 
weapon types throughout the battle. The number of attacking and defending weapon 
types may be different. The following characteristics must be specified for each 
weapon type on each side and do not change over the course of the battle. 

e The number of weapons at the start of the battle. 

e The break point, i.e. the fraction of the starting number of weapons which, if 
reached, would cause the battle to end. For example, if the attacker. would 
withdraw if 50% of a certain weapon type was lost, then the breakpoint for that 
weapon type would be .5. 


e The Lanchester weapon characteristic, 1.e. whether it 1s a square law, linear law, 
logarithmic law weapon, or a hybrid. 


The following characteristics of each weapon type may change from phase to phase of 
the battle. 

e The time required to complete that phase. 

e The rate at which replacements arrive for each weapon type. 


e Attrition coefficients, 1,e. the rate at which a weapon type is attritted by each 
Opposing weapon type in a particular battle phase. 


If there are not more than five weapon types per side, the program prints a dynamic 
display to the M100 screen showing for each weapon type the fraction of starting 
strength that has survived and the breakpoint. Regardless of the number of weapon 
types, the program creates an output file which shows the number of survivors at the 
end of each phase, which weapon reached its breakpoint first, and the number of 


survivors at the end of the battle. 


>The terms attacker and defender are arbitrary and are used only for notational 
purposes. 
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B. REPRESENTING LANCHESTER CHARACTERISTICS OF EACH WEAPON 

The program includes provisions for calculating attrition with traditional 
Lanchester square law and linear law equations or using a Helmbold equation. 
Traditional Lanchester attrition uses different functions for different types of attrition. 
The linear law functions (see Equations 5.1) for changes in strength with respect to 
time are used to model fire that is aimed at a general area in which ‘targets are believed 
to be located An example of a linear law weapon would be artillery fire without 


correction by an observer. 


dx/dt = -axy, and (egress) 
dy/dt = -byx 
The square law functions (see Equations 5.2) are used for fire that is aimed at a point 


instead of a general area. 


dx/dt = -ay, and (egiard:2) 
dy/dt = -bx 
The logarithmic law functions (see Equations 5.3) are used to model non-combat losses 


such as disease. 


dx/dt = -ax, and (eqn 5.3) 

dy/dt = -by 
The reasons for using these equations for modeling these types of attrition are 
explained in [Ref. 9:Chapter 2]. The limitation of using Equations 5.1, 5.2, and 5.3 is 
that they restrict the model to three discrete weapons types. However, there may be 
weapons that do not fall cleanly into any of the three weapon types. For example, 
artillery fire that is corrected by an observer may have a mixture of linear law and 
square law characteristics. The traditional Lanchester equations also assume that the 
full fire power of all the weapons on one side can be brought to bear against all the 
targets on the opposing side. However, in battles where one side vastly outnumbers 
the other, or when there are significant terrain masking or reaction time effects, this 
assumption is not valid. To account for these limitations, R. Helmbold proposed a set 
of modified Lanchester equation in Reference 8. Helmbold’s equations and their 
empirical validation are also discussed in [Ref. 9:pp. 174-181 and Footnote 2.40]. The 
special cases of the Helmbold equations used in this program are those listed as 


Equations 5.4. 
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dx/dt = -a(x/y)?y = -ax2x1-@ (eqn 5.4) 


dy/dt = -b(y/x)®x = -by@x!-® ~ 
They include an additional factor, (x/y) or (y/x)®. If w=1, the result is the 


logarithmic law equation. If @=0, the result is the square law equation. If w=.5, the 
result is an equation which behaves quite like the linear law equation. @ may be set at 
any value between zero and one to characterize weapons which do not fall neatly into 
one of the discrete weapon types. 

In the traditional Lanchester linear law equation dx/dt varies proportionately to 
the number of X survivors and the number of Y survivors. In the Helmbold equation 
when @=.5, dx/dt varies proportionately to the square root of the number of X and Y 
survivors. If dx/dt were equal for the Lanchester linear law and the Helmbold 
equation with @=.5, then aj xy = ary(xy)> Where a; and aj; are the Lanchester and 
Helmbold attrition coefficients respectively. Thus, ayy = ay (xy). Therefore, there 1s 
no single value of ay; that will be equivalent to a value of a; throughout an entire 
simulation because x and y are changing. However, the general shapes of the functions 
Kk, = xy and TET 60) ie are similar enough that they cause attrition to behave about 
the same in both circumstances. A plot of contour lines of dx/dt = 1, 2, and 3 for 
linear law and comparable Helmbold equations is shown in Figure 5.1 where aj;=.02 
and ay =.0002. The contour dx/dt=2 is the same for both formulations.© For 
dx/dt <2 the Helmbold equations produce smaller attrition rates than do the traditional 
Lanchester equations. For dx/dt>2 the Helmbold equations produce larger attrition 
rates than do the traditional lanchester equations. 

The program includes BASIC code for both the traditional Lanchester and 
Helmbold equations. To use the traditional Lanchester equations, lines 223-224 and 
233-234 should be active, and lines 225 and 235 should be commented out or deleted. 
To use the Helmbold equations, lines 225 and 235 must be active, and lines 223-224 
and 233-234 must be commented out or deleted. The weapon type parameters in the 


input file must match the equations which are active in the program. 


©The Helmbold coutours were jittered slightly to avoid being overprinted at 
dx/dt = 2 by the Lanchester contour. 
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Figure 5.1 dx/dt For Linear Law and Helmbold Equations. 


EXPLANATION OF VARIABLES 


¢ AA(NA,ND) is a matrix of the rates at which attacking weapons by type are 


attritted by each type of defending weapon in a particular phase of the battle. 


¢ AB(2,NA) is the matrix for the breakpoints of each attacking weapon type. 


Elements of the first row, AB(1,1), are the breakpoints expressed as fractions of 
the starting total,i.e0 S AB(1,1) S 1 fori=1,2,3...NA. Elements of the second 
row, AB(2,1), are breakpoints expressed as numbers of weapons of type1,1e.0 S 
orem SO) for 1= 1,2,3...NA. 

AT(NA) 1s a vector of tuning parameters that specify the Lanchester 
characteristics of each attacking weapon type. If line 235 1s active and lines 233 
and 234 are commented out, then attrition 1s calculated using the Helmbold 
equation, see equation 5.4. AT(i) is @ and equals 0 for an aimed fire/square law 
weapon, .5 for a weapon similar to a Lanchester area fire/linear law weapon, and 
1 for a source of non-combat/logarithmic law casualties to the defender. If lines 
233 and 234 are active and line 235 is commented out, then attrition 1s calculated 
using standard linear law and square law Lanchester equations, see equations 5.1 


and 5.2. AT(i) equals 1 for linear law weapons and 2 for square law weapons. 
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BB(ND,NA) ts a matrix of the rates at which one item of each attacking weapon 
type attrits each defending weapon type in a particular phase of the battle. 
DB(2,ND) 1s the matrix for the breakpoints of each defending weapon type. 
Elements of the first row, DB(1,j), are the breakpoints expressed as fractions of 
the starting total, 1.e0 & DB(1,)) = Vier, —1,2;3..ND. Blements of the second 
row, DB(2,}), are breakpoints expressed as numbers of weapons of type j, i.e. 0 
= DB(2,}) S SDG) for j= 1,2,3...ND. 

DT(ND) is a vector of tuning parameters that specify the Lanchester 
characteristics of each defending weapon type. If line 225 is active and lines 223 
and 224 are commented out, then attrition is calculated using the Helmbold 
equation, see equation 5.4. DT(j) 1s @ and equals 0 for an aimed fire/square law 
weapon, .5 for a weapon similar to a Lanchester area fire/linear law weapon, and 
1 for a source of non-combat/logarithmic law casualties to the attacker. If lines 
223 and 224 are active and line 225 1s commented out, then attrition 1s calculated 
using standard linear law and square law Lanchester equations, see Equations 5.1 
and 5.2. DT(j) equals 1 for linear law weapons and 2 for square law weapons. 
11,12,f3, and [4 are loop counters. 

MD is the maximum of NA and ND. 

NA ts the number of attacking weapon types. 

ND 1s the number of defending weapon types. 

NI 1s the number of intervals into which a phase of battle is broken. 

NP is the number of phases of battle. 

OA(NA) 1s the vector of horizontal pixel positions for the current screen 
representation of the fraction of surviving attackers by weapon type. 

OD(ND) 1s the vector of horizontal pixel positions for the current screen 
representation of the fraction of surviving defenders by weapon type. 

QA(2,NA) holds the surviving number of each attacking weapon type. QA(1,1) is 
the number at the start of a time increment, DT. QA(2,1) 1s the number after 
attrition by each defending weapon is subtracted. 

QD(ND) holds the surviving number of each defending weapon type after 
attrition by each attacking weapon is subtracted. 

SA(NA) 1s the number of attacking weapons by type at the start of the battle. 
SD(ND) is the number of defending weapons by type at the start of the battle. 
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¢ TF is the termination flag. TF=0 means a breakpoint has not been reached. 
TF=1 means a breakpoint has been reached. 

e TP is the top pixel position for a rectangle or line drawn on the graphical display. 

e TT is the total time that has elapsed in the battle up to the current time 
‘increment. When a breakpoint is reached, or when all phases of the battle are 
completed, TT is time length of the battle reported to LANOUT.DO. 


D. INPUT 
All input to the program is entered into a RAM file, LANIN.DO. The following 
parameters must be entered in the following order and do not change between battle 
phases. 
e NP,NA,ND 
e QA(NA) 
© QD(ND) 
e AB(NA) 
e DB(ND) 
e AT(NA) 
¢ DT(ND) 
The following parameters must be entered in the following order for each phase. 
e TT, NI 
e AR(NA) 
¢ DR(UNR) 
e AA(NA,ND) 
e BB(ND,NA) 
An example of an input file is shown in Figure 5.2 
The situation to be simulated using the parameters in Figure 5.2 is as follows. The 
battle has two phases. 
1. Parameters Common To Both Phases 
The attacker has two weapon types, A,, 1 = 1 and 2. The defender has three 
weapon types, D;, i = 1, 2, and 3. The attacker starts with 200 A,'s and 100 A4’s. 
The defender starts with 100 D1's, 200 D>’'s, and 100 D3’s. The battle will end when 
the first attacking or defending weapon type 1s attritted to 50% of its initial strength, 
1.e. the breakpoint for all weapon types 1s .5. If this simulation is to be run using 
traditional Lanchester equations, both attacking weapons are square law weapons. D, 


is an linear law weapon, and D> and D3 are square law weapons. 


2 


223 
200 100 
100 200 100 


-00012 .014 . 
-00018 .020 . 
-019 .017 
-015 .013 
-O11 .009 

10 50 

so ee 

-4 .% .4 
-00018 .021 . 
.00027 .030 . 
-0285 .0255 
-0225 .0195 
.0175 .00135 





Figure 5.2. Sample Input File. 


2. Situation For Phase One 
The length of the first phase is five hours. The program will break that period 
into 50 segments for computational purposes. Replacements arrive for both attacking 
weapon types at an average rate of one weapon per hour. The replacement rate for all 
three defending weapon types averages .8 weapons per hour. Ay’s are attritted by each 
D, at a rate of .00012 per hour per surviving Al. Aj’s are attritted by each D> at a 
rate of .014 per hour. Rates at which the remaining A; are attritted by each D; are 
listed through the end of the next line. Dy’s are attritted by each A, at a rate of .019 
per hour. Dj’s are attritted by each A> at a rate of .017 per hour. Rates at which the 
remaining D; are attritted by each A; are listed through the end of the next two lines. 
3. Situation For Phase Two 
The length of the second phase is ten hours. The program will break that 
period into 50 segments for computational purposes. Replacements arrive for both 
attacking weapon types at the rate of .5 weapons per hour. The replacement rate for 
all three defending weapon types averages .4 weapons per hour. Ajy’s are attritted by 
each D, at a rate of .00018 per hour per surviving Al. Ay’s are attritted by each D5 


at a rate of .021 per hour. Rates at which the remaining A, are attritted by each D; 


2 


are listed through the end of the next line. Dj's are attritted by each A, at a rate of 
0285 per hour. Dy’s are attritted by each A» at a rate of .0255 per hour. Rates at 
which the remaining D, are attritted by each A; are listed through the end of the next 


J 
two lines. 


EF. OUTPUT 
The program produces two types of output: 

e An output file, LANOUT.DO, which lists the status of each weapon type at the 
end of each phase and at the end of the battle and which weapon types have 
gone below their breakpoints. 

e ea staple sia ay to the screen of the M100 which shows the fraction 
Ofte ashe png ele et ve each micapon type wach has survived until that pice 
The graphical display consists of a rectangle on the M100 screen for each 

attacking and defending weapon type. Each rectangle is 100 pixels wide and five pixels 
high. Each pixel in the horizontal direction represents one percent of the starting 
strength of the weapon type represented by that particular rectangle. In each rectangle 
1s a vertical line showing the breakpoint for that weapon type as a fraction of starting 
strength. As the simulation progresses, another vertical line in each rectangle is 
updated showing the fraction of survivors for that weapon type. The rectangles are 
arranged in two columns, one for attacking and one for defending weapon types. 
Weapon type numbers are printed to the left of the rectangles. If the replacement rate 
drives the number of survivors over the starting strength for some weapon type, the 
vertical line indicating the fraction of survivors will stay at the 100% level and an 
asterisk will be printed to the left of the corresponding rectangle. In addition to the 
rectangles there is a printed line at the bottom of the screen which tells the operator on 


what phase and time interval the simulation is currently working. 


F. EXPLANATION OF PROGRAM COMPONENTS 
A complete listing of this program is at Appendix C. 
1. Initialization Section, Figure 5.3 
Line 120 sets the number of files to two and opens the input file, LANIN.DO. 
Line 121 opens the output file, LANOUT.DO and enters the number of phases in the 
battle, NP, and the number of attacking and defending weapon types, NA and ND. 
Line 122 defines MD as the maximum of NA and ND. If both sides have five or fewer 
weapons types, line 124 sets SF = 1, indicating that the screen of the M100 is big 


enough to handle the graphical display generated during the simulation. If either side 


28: 


"LANCHESTER TIME STEP MODEL 
MAXFILES=2: OPEN'LANIN"'FORINPUTAS1 
OPEN"LANOUT "F OROUTPUTAS2: INPUT#1,NP ,NA,>,ND 


I FNA>NDTHENMD=NAELSEMD=ND 

IFMD<6THENCLS:SF=1 
DIMAA(NA,ND),BB(ND,NA},AT(NA),DT(ND},AR(NA},DROEND } 

DIMQA( 2,NA},QD(ND),AB(2,NA),DB(2,ND},SACNA),SD(ND },OA(NA},OD(ND ) 





Figure 5.3. Initialization Section. 


has more than five weapon types, the graphical display will not be generated. Lines 


130-131 dimension the matrices used in the program. 


2. Entering Common Parameters, Figure 5.4 


‘Enter Initial Quantities of Wpns, Break Points And Wpn Types. 
FORI2=1TONA: INPUT#1,QA(2,I2):SA(I2 J=QA( 2,12): OA(I2 J=127:NEXTI2Z 
FORIZ=1TOND: INPUT#1,QD(I2):SD(I2 J=QD(I2}:0D( 12 J=238:NEXTI2 


FORI2=1TONA: INPUT#1,AB(1,12):AB(2,12 J=AB( 1,12 )¥QA( 2,12 }:NEXTI2 
FORI2=1TOND: INPUT#1,DB(1,12 }:DB(2,12 )=DB( 1,12 )3*QD(I2):NEXTI2 
FORI2=1TONA: INPUT#1,AT(1I2):NEXTI2: FORIZ2=1TOND: INPUT#1,DT(I2):NEXTI2 
TM=0: IFSF=1THENGOSUB6O00 





Figure 5.4 Entering Common Parameters. 


This section enters the rest of the parameters that do not change from phase 
to phase of the battle and initializes the variables controlling the graphical display. 
Line 134 puts the starting quantity of A; into QA(2,1) and SA(i) and sets OA(i) to 127 
which indicates that 100% of A; are surviving at the start of the simulation. Line 135 
performs the same function for D- that line 134 performed for A;. Lines 136 and 137 


J 
enter the fractional breakpoints for A; and D; respectively into AB(1,1) and DB(1,j). 


Lines 136 and 137 also compute the resis in terms of numbers of surviving 
weapons, placing them in AB(2,1) and DB(2,j). Line 138 enters the Lanchester 
characteristic parameters, AT(i) and DT(j), for each weapon type. Line 140 sets TM, 
which keeps track of the time until the end of the battle, to zero. If the graphical 
display to the screen is to be used, line 140 also calls subroutine 600 which sets up the 


Output screen. 
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3. Initialization For Each Phase, Figure 5.5 


FORII=1TONP: PRINT#2,"STARTING PHASE"3I1 
‘Enter Time Spent In Phase Il and # of Intervals 
INPUT#1,TT ,NI:OT=TT/NI 


‘Enter Replacement Rates And Attrition Coefficient Matrices 
FORIZ2=1TONA: INPUT#1 ,AR(I2 ):NEXTI2: FORIZ2=1TOND: INPUT#1 ,DR( I2):NEXTI2 
FORI2=1TONA: FORI3=1TOND: INPUT#1 ,AA(I2,13 ):NEXTI3:NEXTI2 
FORI2Z=1TOND: FORI3=1TONA: INPUT#1,BB(I2,I3):NEXTI3:NEXTI2 





Figure 5.5 Initialization For Each Phase. 


Line 143 sets I1, the phase counter, and prints a heading to the output file. 
Line 146 enters the length of the current phase, TT, and the number of intervals, NI, 
into which the phase will be broken and computes the length of each interval, DT. 
The choice of NI is a compromise between two competing objectives: accuracy and 
time required to complete the simulation. As NI becomes larger, DT becomes smaller 
and the simulation more closely approximates the continuous, mutual attrition that is 
the basis of Lanchester equation theory. If NI is small and DT is large, then the 
simulation tends to discount the attrition which takes place during a phase because the 
program assumes force levels are constant throughout an interval, DT. However, if NI 
is to large, the time required to run the simulation increases linearly. The relationship 
between accuracy and speed is also a function of the attrition coefficients and number 
of weapon types on each side. 

Line 152 enters the replacement rates for each weapon type. Lines 154-156 
enter the attrition coefficients for each weapon type. 

4. Attrition Calculations For A Phase, Figure 5.6 

This section breaks a phase into NI intervals of length DT, calculates the 
attrition during each interval, and tests whether that attrition has caused some weapon 
type to reach its breakpoint. 

Line 202 sets the interval counter, 12, adds the length of the interval to the 
length of the battle, TM, and prints a message to the screen telling the operator what 
phase and interval is currently being processed. 

Lines 222-227 calculate the attrition to attackers based upon the quantities of 


each weapon surviving at the start of the interval. QA(2,1) holds the current quantity 
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‘Fight Phase Il. 

FORI2Z2=1TONI: TM=TM+DT:PRINT9241,"Phase:"3I13", Increment"3I2Z3;"out of"3NI 
‘Fight Time Increment DT. 

‘Update number of attackers 

FORI3=1TONA: QA(1,I3 J)=QA( 2,13 ):NEXTI3: FORIS=1TONA: FORI4=1TOND 

IFDT(IG¢ J=LTHENQA( 2,13 J=QA( 2,13 )-AA(I3 >I4 )¥QD( 14 )¥QA( 2,13 )¥DT:GOTO2Z26 
QA( 2,13 J=QA( 2,13 )-AACI3 »I4¢ )¥QD( IG )*DT 

"QA( 2,13 J=QA( 2,13 )-AA( 13,14 }¥(QA( 2,13 )/QD( 14 ) J“ DT( 14 )*¥QDI 14 )¥dt 
NEXTI4: QA( 2,13 J=QA( 2,13 )+AR(I3 J¥DT: IFSF=1THENGOSUB650 

NEXTI3 

‘Update number of defenders 

FORI3=1TOND: FORIG=1TONA 

IFAT(IG¢ )=1THENQD( 13 J=QD(1I3 )-BB( 13,14 )*QD(I3 )¥QA( 1,14 )*DT: GOTO236 
QD(I3 J=QD(I3 )-BB( 13,14 )¥QA(1,I4¢ )¥DT 

'QD( 13 J=QD( 13 )-BB( 13 ,14)*( QD(I3 J/QA( 1,14) JA ATIIG )¥QA( 1,14 )¥dt 
NEXTI4¢: QD(1I3 J=QD(I3 )+DR(1I3 )*DT: IFSF=1LTHENGOSUB660 

NEXTI3 

GOSUB300 : NEXTI2 

IFI1L=NPTHENGOSUB350:CLS:PRINT"Output is in file LANOUT.DO.":END 
PRINT#2,"Status After Phase"3I1:GOSUB361:NEXTI1 





Figure 5.6 Attrition Calculations for a Phase. 


of A; and is therefore decremented by lines 222-227. Attrition to each D; should, for 
consistency with the attrition to the A,’s, also be based upon the quantity of A,’s 
surviving at the beginning of the interval. Therefore, the quantity of each A; surviving 
at the beginning of the interval is saved by line 222 in QA(1I,i) for use during the 
attrition calculations for D,. Line 222 also sets the A; counter, I3, and the D; counter, 
14. 

If the simulation is to be done with traditional Lanchester linear law and 
square law equations, lines 223-224 must be active and line 225 must be commented 
out or deleted. If the simulation is to be done with a Helmbold equation, lines 223-224 
must be commented out or deleted and line 225 must be active. The program 
displayed in Figure 5.6 has the Helmbold equations commented out. Line 222 
calculates the attrition of A; by D; based upon a linear law, Equation 5.1. Line 223 
calculates attrition based upon the square law, Equation 5.2. Whether the linear law 
or square law is used is based upon AT(i) and is determined in line 222. If line 225 
were active, it would calculate attrition using the Helmbold equation, Equation 5.4. 
After attrition by each D; is calculated, line 226 adds the replacements for A; received 
during the interval. Line 226 also calls subroutine 650 which updates the graphical 


display to reflect the attrition to each A.. 
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Lines 232-237 calculate the attrition to defenders based upon the quantities of 
each weapon surviving at the start of the interval. QD(i) holds the current quantity of 
D;. Line 232 sets the D; counter, 13, and the A; counter, 14. If the simulation is to be 
done with traditional Lanchester linear law and square law equations, lines 233-234 
must be active and line 235 must be commented out or deleted. If the simulation is to 
be done with a Helmbold equation, lines 233-234 must be commented out or deleted 
and line 235 must be active. The program displayed in Figure 5.6 has the Helmbold 
equations commented out. Line 232 calculates the attrition of D; by A; based upon a 
linear law, Equation 5.1. Line 233 calculates attrition based upon the square law, 
Equation 5.2. Whether the linear law or square law is used is based upon DT(i) and is 
determined in line 232. If line 235 were active, it would calculate attrition using the 
Helmbold equation, Equation 5.4. After attrition by each A; 1s calculated, line 236 
adds the replacements for D; received during the interval. Line 236 also calls 
subroutine 660 which updates the graphical display to reflect the attrition to each D,. 
Line 240 calls subroutine 300 to check for whether a breakpoint was reached during 
the interval. If no breakpoint was reached, a new interval is begun. 

If all the intervals in the last phase are completed without reaching a 
breakpoint, line 242 calls subroutine 350 which prints the status at the end of the battle 
to the output file. If the current phase is not the last phase, then line 245 prints a 
header to the output file, calls subroutine 361 which prints the status at the end of the 
current phase to the output file, and starts the next phase. 

5. Breakpoint Subroutine, Figure 5.7 

This section determines whether any weapons have reached their breakpoints, 
and, if so, prints that information in the output file, and ends the program. Line 320 
sets TF = 0, indicating that no breakpoints have been reached. Line 320 then starts a 
loop which tests whether any A; have reached their breakpoints. If so, then lines 
322-324 print a message to the output file specifying the weapon type, the breakpoint, 
and the quantity of that weapon type that survived. Lines 335-340 test whether any D; 
have reached their breakpoints, and if so, print a message to that effect to the output 
file. If no breakpoints have been reached, then line 340 returns control to the main 
program to begin attrition calculations in the next time interval. 

The default battle termination criterion is that at least one weapon type must 
be below its individual breakpoint at the end of a time interval. However, the operator 


may wish to edit the program before running it, adding more sophisticated termination 
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‘Check Whether Breakpoint is reached. 

TF=0: FORI3Z=1TONA: IFQA( 2,13 )>AB( 2,13 JTHEN3S25 

TF=1: PRINT#2,"Attacker Wpn"31I33;"Is Below Breakpoint" 
PRINT#2," Bp ="3:PRINT#2,USING"BHH#.HH"3AB( 2,13 )3 

PRINT#2," Current Level ="'3: PRINT#2,USING"H##H. #H#"3QA( 2,13) 
NEXTIZ 

FORI3=1TOND: IFQD(1I3 J>DB( 2,13 JTHEN340 

TF=1:PRINT#2,"Defender Wpn"3;I3;"Is Below Breakpoint" 
PRINT#2," Bp ="3:PRINT#2 ,USING' HHH. #4" 53D0B( 2,13 )3 

PRINT#2," Current Level ="3:PRINT#2 USING" H###H. #4" 5 QD( 13 ) 
NEXTI3: IFTF=OTHENRETURN 

PRINT#2,""': PRINT#2,"":PRINT#2,"SUMMARY AT END OF BATTLE" 
PRINT#2,"":PRINT#2,"Time Elapsed During Battle ="; 

PRINT#2 ,»USING'H#RF. #24" 37TM: PRINTHZ,"": GOSUBS61 ° 
CLS:PRINT"Output is in file LANOUT.DO":END 

PRINT#2," Att Hpn Breakpoint Current Level" 
FORI3Z=1TONA: PRINT#2 ,»USING'HHHRHH" 5133 . 

PRINT#2 ,USING"##HHHHEHHH . HH" 5 AB( 2,13 )3;QA( 2,13 ):NEXTI3:PRINT#2,"" 
PRINT#2,"" Def Hpn Breakpoint Current Level" 
FORI3=1TOND: PRINT#2 ,USING'H#H#HH" 5135 

PRINT#2 ,USING' HRERHHHHEH . HH" 5 DBI 2,13 )3QD( 13 ): NEXTI3:PRINT#2,""; RETURN 





Figure 5.7 Subroutine To Test For Breakpoints. 


criteria to the default criteria. For example, if the operator wants the battle to 
terminate when A, or A, reach half their starting strength or when A, reaches 60% 
and A» reaches 70% of their starting strength the operator should: 

e Put .5 as the individual breakpoints for A, and A, in the input file and 


e Put the program lines shown in Figure 5.8 into the program after line 325. 


330 IFQA(2,1)/SAC1)>. 60RQA(2 52 D/SA 2 >. 7THEN335 


331 TF=1:PRINT#2,"Special Termination Criterion Met" 





Figure 5.8 Example Of An Additional Termination Criterion. 


If a breakpoint has been reached, then lines 350-368 print the status of both 
sides at the end of the battle. That end of battle status report includes a header, time 
elapsed during the battle, and a list of attacking and defending weapon types with their 
breakpoints and number of survivors. Lines 361-368 are called as a subroutine from 
line 352 because lines 361-368 are also used to print the summary at the end of each 


phase and can therefore not terminate the program. 
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6. Graphics Display Initialization Subroutine, Figure 5.9 


'Set up output screen : 

PRINT"Wpn # Attacker Defender" 
FORI1=1TOMD: PRINTUSING'##" 511 

TP=2+1I1*8 

IFI1>NATHEN630 


LINE(18,TP J-(119,TP+4),1,8:BP=18+INT(100*AB(1,I1)) 
LINE(8P-1,TP+1)-(BP,TP+3),1,8 

IFIL>NDTHEN635 

LINE(138,TP J-(239,7TP+4),1,8:BP=138+INT( 100*DB(1,1I1)) 
LINE(BP-1,TP+1)-(BP,TP+3),1,8 

NEXTI1: RETURN 





Figure 5.9 Graphics Display Initialization Subroutine. 


The graphics display consists of a rectangle on the M100 screen for each 
attacking and defending weapon type. Each rectangle is 100 pixels wide and five pixels 
high. Each pixel in the horizontal direction represents one percent of the starting 
strength of the weapon type represented by a particular rectangle. The rectangles are 
arranged in two columns, one for attacking and one for defending weapon types. This 
subroutine draws the rectangles, labels the columns “Attacker” or “Defender”, puts a 
vertical line in each box at the breakpoint for that weapon type, and labels the rows of 
boxes with the weapon type number. | 

Line 610 prints the header on line one of the screen. Line 620 starts a loop 
that writes the weapon type number, I1, and prints the rectangles; Line 623 calculates 
the vertical pixel position, TP, for the top of the rectangles of weapon type Il. Line 
625 tests whether to draw a rectangle next to Weapon type number I] in the “Attacker” 
column. If so, the first statement on line 627 draws the rectangle. The second 
statement on line 627 calculates the horizontal pixel location in the rectangle of the 
breakpoint for that weapon type. Line 628 draws a double line at the breakpoint in 
the rectangle. Lines 630-635 test whether a rectangle should be drawn in the defender 
column. If so, they draw the rectangle and insert the breakpoint in the same manner 
as was done in lines 620-628. Line 635 returns control to the main program when 


there are no more rectangles to be drawn. 
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7. Updating The Graphical Display, Figure 5.10 


"Update screen output of attackers 

TP=3+13*8 

LINE(OA(I3),TP )-(OA(I3),TP+2),0 

OA(I3 J=18+INT(100*QA( 2,13 )/SA(I3)) 

IFOA(I3 )>LISTHENOA(I3 J)=118:PRINTO(13*40+2),"x"': GOTO659 
PRINTOI3*40+2," " 
LINE(OA(I3),TP)-(OA(I3),TP+3),1:RETURN 


‘update screen output of defenders 

TP=3+1I3%*8 

LINE(OD(I3}),TP)-(OD(I3),TP+2),0 

OD( 13 J=138+INT(100*QD(I3 )/SD(1I3)) 

IFOD( 13 )>238THENOA( 13 )=238: PRINTOI3*40+22 »"*" 3: GOTO669 
PRINTOI3*40+22," " 

LINE(OD(I3),TP )-(OD(I3),TP+3),1:RETURN. 





Figure 5.10 Subroutines To Update The Graphical Display. 


This section includes two subroutines which update the vertical line in each 
rectangle which indicates the fraction of survivors for that weapon type. The 
subroutine in lines 650-659 updates the attacker rectangles; lines 660-669 perform the 
same function for defender rectangles. Line 653 sets TP, the location of the top pixel 
of the vertical line for Ay3. Line 655 erases the old vertical line, the horizontal pixel 
position for which was stored in OA(I3). Line 656 calculates the horizontal pixel 
position for the new vertical line based upon the fraction of the starting strength of 
Ay3's which currently survives. If the reinforcement rate exceeds the attrition rate and 
drives the number of survivors over the starting strength for Aq3, line 657 holds 
horizontal position of the vertical line at the 100% level and prints an asterisk next to 
the corresponding rectangle. If the number of survivors is less than the starting 
strength, line 658 prints a blank space next to the corresponding rectangle. Line 659 
writes the new vertical line to the screen showing the fraction of Ay;3’s which survive. 


Lines 660-669 perform the same function for D; that lines 650-659 perform for A;. 


G. EXAMPLE SIMULATIONS 
1. Example #1 
The first example uses the input file shown in Figure 5.2 The output file for 


that simulation is shown in Figure 5.11. 
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STARTING PHASE 1 
Status After Phase 1] 


Att Won Breakpoint Current 
1 100.00 17s. 
2 50.00 68. 
Def Wen Breakpoint Current 
1 50.00 79. 
2 100.00 184. 
ES 50.00 89. 


STARTING PHASE 2 


Attacker Wpn 2 Is Below Breakpoint 
Bp = 50.00 Current Level = 48.85 


SUMMARY AT END OF BATTLE 

Time Elapsed During Battle = 7.20 

Att Wpn Breakpoint Current 

1 100.00 157. 

2 50.00 48. 

Def Wen Breakpoint Current 

50.00 66. 

100.00 174. 

50.00 84. 





ricwune ole Outour Piles LANOUT.DO, ForExampler7 |: 


2. Comparing The Lanchester and Helmbold Linear Law Equations 
Examples 2 and 3 compare the differences between using a traditional 
Lanchester linear law equation (Example 2) and a Helmbold equation (Example 3). 
The scenereos for Examples 2 and 3 share the following elements. 


e The battle has only one phase with a maximum length of 10 which 1s broken into 
100 intervals. 


° ae sides have two weapon types. Each weapon type has a starting strength of 


e The breakpoints for all weapon types are 50%. 


e All weapon types are linear law weapons. Attrition is_calculated using Equation 
5.1 for Example 2 and using Equation 5.4 (@=.5) for Example 3. 


e There are no replacements for any weapon type. 
The only differences in the scenereos for Examples 2 and 3 are the attrition 
coefficients. 
a. Example #2 
The input and output files, LANIN.DO and LANOUT.DO, for Example 2 
are in Figure 5.12. Since these attrition rates are for the Lanchester linear law, the 


dimensionality of the rates for the attacker are (number of attacker casualties) per 
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(number of attackers) per (number of defenders) per (unit time). The dimensionality of 


the rates for the defender are the same with the rolls reversed. 


Input File: Output File: 


ee 2 STARTING PHASE 1 
100 100 Attacker Wpn 1 Is Below Breakpoint 
Bp = 50.00 Current Level = 49.40 


SUMMARY AT END OF BATTLE 
Time Elapsed During Battle = 5.20 


Att Wpn Breakpoint Current Level 
z ‘ 49.40 
2 .50. 62.53 
-00075 .00075 
-0005 .0005 ; Def Wen Breakpoint Current Level 
-00025 .00025 i 50.00 82.17 
-00025 .00025 2 50.00 82.17 





Figure 5.12 Input And Output Files For Example #2. 


b. Example #3 

The input and output files; LANIN.DO and LANOUT.DO, for Example 3 
are in Figure 5.13. The attrition rates in this example are for the Helmbold Equation 
with @ = .5, the Helmbold equivalent of the Lanchester linear dimensionality of the 
rates for the attacker are (number of attacker casualties) per (attacker)> per 
(defender)*» per (unit time). The dimensionality of the rates for the defender are the 
same with the rolls reversed. 

To generate Helmbold coefficients that are comparable to the Lanchester 
linear law coefficients, the Lanchester coefficients must be adjusted by the difference in 
the dimensionality, 1.e. multiplied by [(number of attackers)(number of defenders)J°>. 
In this example it means multiplying the Lanchester coefficients by 100. 

The results of Examples 2 and 3 show good agreement. 

e The simulation using the Helmbold equations ended about 21% faster than did 
the battle using Lanchester equations. The same weapon type reached its 
breakpoint first in both cases. 


e The BLES between the number of survivors for other weapon types was 
quite small. 
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Input File: Output File: 


STARTING PHASE 1 
Attacker Wpn 1 Is Below Breakpoint 
Be = 50.00 Current Level = 49.85 


SUMMARY AT END OF BATTLE 
Time Elapsed During Battle = 4.10 


Att Wpn Breakpoint Current Level 
1 50.00 49.85 
2 50.00 64.67 
-0000075 .0000075 
-000005 .000005 Def Wen  Breakpoint Current Level 
-0000025 .0000025 1 50.00 82.79 
-0000025 .0000025 2 50.00 82.79 





Figure 5.13. Input And Output Files For Example #3. 
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VI. GEQMETRIC PROGRAMMING 


A. GENERAL 
The program described in this chapter solves nonlinear programming problems in 
which: 
e The objective function is to be minimized. 
e The objective function and constraints are posynomials. 
e The number of terms,’ T, minus the number of variables, N, must equal one. 


© The coefficients,® c of all terms must be strictly positive. 


m,t? 
e All components of the vector of decision variables, ¥, must be strictly positive at 
optimality. 
e Constraints must have the form of a posynomial on the left hand side that is less 
than or equal to one. 
Geometric programming has the distinctive feature of calculating the optimal 
value of the objective function before the optimal values of the decision variables are 
calculated. Geometric programming also produces weights, O1, t= 1,2,3,...,1T, associated 


with each term. For example, in applications where the c, are prices and the objective 


’The number of terms in the objective function and in the constraints. 


ST wo subscripting systems are used throughout this chapter. The first uses the 


letters mand t where ¢ = 1,2,3,...,T,, 1s the number of the term in the mth posynomial. 


m 
m=0 refers to the objective function; m= 1,2,... refers to the constraint numbers. T 


m 
is the number of terms in the mth constraint. When problems of only one posynomial 
are being discussed, the m is omitted. The first system also includes the letter n to 
identify components of the decision variable vector, ¥, where n= 1,2,...,.N. The second 
system numbers the terms without starting again at 1 at the beginning of each 
constraint. Each term is numbered t’, t’=1,2,...,T’ where T’ is the number of terms in 
the objective function and the constraints. The first term of the objective function is 
denoted by t’ =1. The other terms in the objective function are then numbered from 
left to rnght. Then the terms in each constraint in turn are numbered from left to nght. 
For example, in Figure 6.1, t=1, m=0 (or t’=1) refers to 40x,x5 and t=2, m=1 (or 


t = 5) refers to (Oe oe ee 
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function minimizes total cost, the 6, for objective function terms are the proportion of 
cost that term t contributes to optimal total cost, fy). These weights are invariant 
with respect to the prices, c,, associated with each term. 

1. Definition Of A Posynomial Function 


The function f(x) is posynomial if it has the form 


T N q 
a enP, co) where p(y) =e Ot (eqn 6.1) 
t=] n=l 


where 
e JT is the number of terms. 
* c, are positive scalar constants. 


° x is the vector of decision variables, (X1,X>,.--:Xx)- 


e The only restriction on the exponents, a is that they be real numbers. 


m,n,t? 
A posynomial differs from a polynomial in that the coefficients of a posynomial must 


be strictly positive and its exponents, a,, ,, need not be positive integers. 


n,v’ 
An example of a problem meeting these conditions is in Figure 6.1. 


Min 40x4X9 7 20x9Xx3 
Subject to: 


, Pig ee ~ iy ea al 
mee 0, 1=1,2)3 





Figure 6.1 Geometric Programming Problem In Standard Form. 


Geometric programming solves a problem of this kind by solving its dual. 
When, as specified above, the number of terms minus the number of variables equals 
one, then the problem has a unique solution. T - N - 1 is by convention called the 
degree of difficulty. If the degree of difficulty 1s greater than zero, then another 
nonlinear program must be solved to find the optimal 5, While this new nonlinear 
program is frequently easier to solve than the original problem, its solution 1s beyond 
the scope of this chapter which is limited to problems with a degree of difficulty of 


ZeCTO. 
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B. MATHEMATICAL BASIS FOR GEOMETRIC PROGRAMMING 
The mathematical basis for geometric programming is summarized in [Ref. 10:pp. 
494-522]. and explained in detail in Reference 11. The following explanation is 
provided for tutorial purposes and is an adaptation of the explanation in [Ref. 10:pp. 
496-502]. Notation in this chapter is consistent with that used in Reference 10. 
a. Unconstrained Minimization Of Posynomial Functions 
Historically, geometric programming has been based upon and took its 


name from the arithmetic-geometric mean inequality: 


T oes T 

y On lly, ‘ if v,,0, = Omand y 0 eae (eqn 6.2) 

t=1 t=1 t=1_ 
The equality holds only when v;=v5=....,vy. If u, is defined as o —8,.e5 then 
Equation 6.2 becomes 

T T 5 

yu, = Wie (eqn 6.3) 

t=1 t=1 


The equality holds if 6, = te U,- Let u, be a posynomial term as described in 
Equation 6.4. 


a a 
We = ¢ be OoinS cae x n,t (eqn 6.4) 
A posynomial function, f{¥), is given by Equation 6.5. 
T 
f(x) = vu, (eqn 6.5) 
t=1 
When the posynomial terms are substituted into Equation 6.3, the inequality becomes 
T T N an 5, 
»» u, = I {Cc, Wx, © 1/6.) (eqn 6.6) 
t=1 t=1 n=1 
or 


where @ is the sum over t of an tt: 


T Oo. N 
u, 2 {IT Ce,/5,1 {11 x,%} (eqn 6.7) 
1 t=1 n=1 
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Since the only restriction on 6, has been that they be positive and sum to 
one, they may be chosen such that » = 0 for n = 1,2,...,N. If this selection is made, 


Equation 6.5 becomes 
>» eype(t) 2 n 1 (e/8)"" (eqn 6.8) 


Since equality holds when 6, = u,/})u,, then 


T iy 5 
min > u, = max I (c,/6,) : (eqn 6.9) 
t=1 t=1 


D8, = 1 and Yan’ , = 0 forn = 1,2,...,N 


Therefore, the minimization of the posynomial function is the same as the 
maximization of the nonlinear function in Equation 6. subject to linear constraints. 
The linearity of the constraints means that 6,” * and va can be computed with linear 
algebra as explained below, instead of with nonlinear programming. These 
minimization and maximization problems are duals. Since equality holds if and only if 


6 


p= i ae it follows that the relationship between optimal values of 6 and x is 


or) 


= {px AQ) oF (eqn 6.10) 


N *% *% 
6,” = fc, Mx, } Bmx’). (eqn 6.11) 
n=l 


If the degree of difficulty is zero, then the matrix of exponents, an with 
another row of |’s appended to the top makes a square matrix. Rows of the exponent 
matrix correspond to variables and columns correspond to terms. The row of I's 
corresponds to the constraint that the sum over t of 6, equals 1. The 5” can be 
obtained by solving a set of T simultaneous linear equations A& = b. The first 
element of the b vector is | and the remaining elements are 0. The optimal value of 
the objective function, yx") can then be obtained by inserting 6. into Equation 6.12. 
Equation 6.12 is based upon Equation 6.9. 


Tt 


# «Oy 
i) = Sea ) (eqn 6.12) 
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Finally, the optimal values of the decision variables, x, can be determined by solving 


nN?’ 
the set of T equations of the form specified in Equation 6.13. 


N P 
Ds an tin(xy ) = Int y )8, jc.) fort = haere. (eqn 6.13) 
n=1 


This set of equations is overconstrained since there are only T-1 decision variables. 
Therefore, only T-1 equations at required. Solving these T-1l equations 
simultaneously produces p, = ag - which are then converted to the Dat values 
of the decision variables by x, = e ™. 
b. Inequality Constraints 

This section discusses the addition of posynomial inequality constraints to 
the unconstrained problem discussed above. For notational purposes the objective 
function and constraints will be numbered m=0,]1,2,3,....M. The objective function 1s 
designated m=0, and the constraints are designated m=0,]1,2,3,....M. A primal 


constrained posynomial would have the form 


T N ay 
Min > {cy, IT x, st) (eqn 6.14) 
t=1 | n=l 
Sulpject to, 
T N 4 
es) = = orn Bees mnt < |] form = 1,2,...,.M (eqn 6.15) 
=1 n=1 


where xX, > Omen = 142.5.44No do t are the weights for the terms in the objective 


function, then 


*k *k %S 
80 t es Leo Po (X )41/fo(x ) for t= 1,2,3,..., 1g. (eqn 6.16) 


If hon are the Lagrange multipliers associated with constraint m, then 


S 
hon = 2 Orm,t and (eqn 6.17) 
om tm = Cnt Pm, t(%)- (eqn 6.18) 
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The dual geometric program is 


MOT 
r 5 
Max il Ii Leth! Om.t2 ray (eqn 6.19) 
m=1 t=1 | 
Subject To: 
T 
» 84 = | (eqn 6.20) 
t=1 
fe 
> > om: — Ofori = 12)3,,N. (eqn 6.21) 
m=1 t=1 
Ue 
| ns \? omit (eqn 6.22) 
t=1 
and ont hon ae 


The 6,,; are calculated using Equations 6.20 and 6.21 as a set of 
simultaneous linear equations. The optimal value of the objective function 1s 
calculated by multiplying the unconstrained optimum by IN(A,) mas in Equation 
6:23: 


* i 6, , hin 
fold) = TE (cp/5,) ° TH Om) (eqn 6.23) 
=1 m=1 
4 is calculated using Equations 6.24 and 6.25. 
a % % - 
Y 29 nln(Xq ) = In(fg(x")89.¢ /o,:) fort = 1,2,..Tg, and (eqn 6.24) 
n=1 
N * *% * 
rN mel Omn . Ent Am ) (eqn 6.25) 
n=1 


Pome 2,-7.2. 3M, and Se 0. 2 
As with the unconstrained problem these equations are linear in tao) =f) oeiler 
solving for p,, using Equations 6.24 and 6.25 as a set of simultaneous linear equations, 


se : p 
the decision variables are calculated using x, = e as 
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D. 


EXPLANATION OF VARIABLES 
BIC(NT + 1,NT*2) holds the A matrix in the subroutine which solves simultaneous 
linear equations of the form Ax=b. 

B2(NT) holds the b vector in the subroutine which solves simultaneous linear 
equations of the form Ax=b. 

B3(NT,NT-1) stores the exponents of the variables in each term. Each row of 
the matrix corresponds to a term; each column corresponds to a variable. Ifa 
variable is not stated explicitly in a term, then its entry in this matrix is zero. 
CT(3,NC,MN) holds three values for each term. CT(1,m,t) holds the coefficient, 
Cnt CT(2,m,t) holds the weight, Ome 


for each term. m=O refers to terms in the objective function. m=1,2,....NC 


for each term. CT(3,m,t) holds p,, () 


refers to terms in the mth constraint. t= 1,2,...,.NT(m) specifies a particular term 
in the objective function or a constraint. 

FS is the optimal value of the objective function. 

11,12, and [3 are loop counters. 

K2,K3,K4,....K9 are variables in the simultaneous linear equation solving 
subroutine. This subroutine is documented in Appendix E. 

LM(NC) holds values of han m= 1)2,--.,.NC where a misathessumauel Ome for the 
mth constraint. Ag = I. 

MN 1s the maximum number of terms in any constraint or the objective function. 
NC is the number of constraints. 

NT(NC) is the number of terms in each constraint. 

NT is the number of terms in the objective function and the constraints. 

NV is the number of decision variables, i.e. the number of components in the 


VEClOr ¥. 


INPUT 
Problem parameters are entered into an input file, GEOIN.DO, before the 


program is executed. GEOIN.DO must contain the following parameters in the order 


specified. 


The number of terms, NT, the number of variables, NV, and the number of 
constraints, NC. 

The number of terms in the objective function, NT(0), and the number of terms 
in each constraint, Ni(m) m= lize, 


The coefficients of each term, ¢ CY(1,m,t) m=0,1,2,2.,NGa— 2 Na 


m,t’ 
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e The matrix of exponents, a for each variable in each term. Row n of the 


n,t? 


matrix corresponds to the variable x, n=1,2,...,N. Column t of the matrix 


n, 
corresponds to the term p,{X), t'=1,2,...,T’. 


The input file for the problem specified in Figure 6.1 is in Figure 6.2. 


Problem: 
Min 40x1X> + 20x5x3 


Subject To: 


2X, 


——*, 
X- 0 





Figure 6.2 Sample Input File, GEOIN.DO. 


E. OUTPUT 
The program prints the following output to the screen. 
¢ The optimal dual variables, ne 
¢ The optimal value of the objective function, Ry). 
e The value of each Baa (OMe 


¢ The optimal value of each component of x, i 


F. EXPLANATION OF PROGRAM COMPONENTS 

A complete program listing is located at Appendix D. 

1. Initialization And Input, Figure 6.3 

Line 110 opens the input file, GEOIN.DO. Line 120 enters the number of 

terms, NT, and the number of variables, NV, from GEOIN.DO, and sets K9 equal to 
NT for use in the simultaneous linear equation solving subroutine. Line 122 checks 
whether the degree of difficulty is equal to zero. If it is not, an error message is printed 
and the program ends. Line 130 enters the number of constraints, NC, and dimensions 
the vector NT(NC), which holds the number of terms in the objective function and in 
each constraint. Line 140 enters NT(m), m=0,1,2,....NC, and computes MN, the 
maximum over m of NT(m). Line 143 dimensions the matrices required for the 


program. Line 145 enters the coefficients c placing them in CT(1,m,t). 


m,t? 
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‘Geometric Programming Program 

OPEN'GEOIN"FORINPUTAS1 

INPUT#1 NT »NV: K9=NT 

IFNT-NV<>1LTHENPRINT"***XERROR: Degree of Difficulty <> 0":END 
INPUT#1 NC: DIMNT (NC ) 

MN=0: FORI1=OTONC: INPUT#1,NT(I1):IFNT( Il J>MNTHENMN=NT(I1):NEXTI1 


DIMCT(3»NC »MN) »>LMENC ),BLENT4+1,NT¥2),B20NT) > BSENT NV) 
FORI1=OTONC: FORI2=1TONT(I1): INPUT#1,CT(1,11,12):NEXTI2:NEXTI1 
FORT1=1TONT(O):B1(1,I21)=1:NEXTI1 

FORII=NT(0 )4+1TONT:B1(1,11)=O0:NEXTI1 

FORI1=2TONT : FORI2=1TONT : INPUT#1>,B1(11,12):B3(1I2,11-1)=B1l(1I1,1I2) 
NEXTI2:NEXTI1 





Figure 6.3 Initialization and Input. 


Lines 150-162 enter elements of the A matrix required to solve the 
simultaneous linear equations Ad=b into B1(,). Equations 6.20 and 6.21 are the basis 
for this set of simultaneous linear equations. Lines 150-155 fill the first row of BI with 
Ones in columns corresponding to objective function terms and zeros in columns 
corresponding to other terms. B1(1,t) corresponds to Equation 6.20. Lines 160-162 
enter the exponents of variables in each term into Bi and B3. In B1 rows 2 through 


NV+1=NT correspond to variables x, and columns | through NT correspond to 


n 
terms t’=1,2,...,T’. Storing the exponents in B3 is necessary because the simultaneous 
linear equation subroutine changes the matrix in Bl, and the exponents are required 
for later calculations. B3 is the transpose of Bl because of the nature of the 
calculation for which B3 is later recalled. 


2. Calculating The Weights For Each Term, Figure 6.4 


PRINT": PRINT'*X*COMPUTING DELTA'S*x" 
B2(1 )=1: FORI1=2TONT :B2( 11 )=O:NEXTI1 
GOSUB 9800 


CLS: I1=1:FORI2=0TONC: FORIS=1TONT(I2):CT(2,12,13)=B1l(1I1,1) 
PRINT"DELTAC™3I23"5"%5135") = "3 
PRINTUSINGHHRHF .RRRE' SCT 2,125,135): L1L=L1+1:IFLI1>STHENGOSUB6OO 
NEXTI3:NEXTI2:GOSUB600:CLS 





Figure 6.4 Calculating 6_, ,. 


72 


Line 172 places the simultaneous linear equation b vector into B2. B2(1)=1 is 
the right hand side of the constraint in Equation 6.20. The right hand sides of the 
constraints based upon Equation 6.21 are zero. Line 180 calls the simultaneous linear 
equation solver in subroutine 9800 which leaves 5° in Bic Yt — eee ere ines 
200-205 place 5” into CT(2,m,t) and prints 0, i to the screen. 

3. The Optimal Objective Function eine, Fane 6.5 


PRINT": PRINT" COMPUTING OPT OBJ FN VALUE*x" 
FORI1=OTONC: LM( 11 )=0: FORIZ=1TONT(I1):LM(I1 J=LM(I1)+CT(2,11,12) 
NEXTI2:NEXTI1 


FS=1: FORI1=OTONC : 
FORI2Z=1TONT(1I1):FS=FS*(CT(1,I1,I2)/CT(2,I1,12) )ACT(2,11,12) 
NEXTI2:FS=FS*(LM(I1)“LM(I1)):NEXTI1 ; 

PRINT": PRINT" FX =''3 ;PRINTUSING'#### .#4#" 5 FS :GOSUB600:CLS 





Figure 6.5 The Optimal Objective Function Value. 


Lines 212-214 compute A, which equal the sum over t of om t for every 
constraint and the objective function. Lines 220-224 compute the optimal objective 
function value based upon Equation 6.23. Line 229 prints the optimal objective 
function value. 


4. Optimal Decision Variable Values, Figure 6.6 


‘Compute optimal x(n) 

K9=K9-1 

FORI1=1TOK9: FORI2=1TOK9:B1(1I1,I2 J=B3(I1,I2):NEXTI2Z:NEXTI1 
CC=1:FORI1=OTONC: FORIZ=1TONT(I1) 

CT(3,I1,I2 J=(CT(2,I1,I2)/CT(1,I1,12)/LM(I1)) 
IFII=OTHENCT(3,I1,12J=CT(3,I1,I2 )*FS 

B2(CC J=LOG(CT(3,I1,12)):CC=CC#1:NEXTIZ:NEXTI1 

PRINT" P(m,t)* = opt. value of term t, constr. m» divided by its coef." 
FORII=OTONC: FORIZ=1TONT(I1):PRINT"P(''5I1s","sIZ23" J" =" 5 
PRINTUSING#### . ####"5CT(3,I1,I2):NEXTI2:GOSUB600:NEXTI1:CLS 
PRINT" :PRINT"** Computing Opt Values Of X(n) *%" 

GOSUB 9800 

CLS: FORI1=1TOK9:PRINT'X*("3I13") = "3 
PRINTUSING"###8884 . FHHHH" SEXP(B1I(I1,1)) 

IFIL>STHENGOSUB600 

NEXTI1:GOSUB600: END 





Figure 6.6 Optimal Decision Variable Values. 
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After 6 and the optimal value of the objective function, fox), have been 
calculated, this section calculates y . The basis for these calculations is Equations 6.24 
and 6.25. The section solves a set of simultaneous linear equations A Cin(y )J =b and 


then solves for y . 


Since the number of variables is one less than the number of terms, line 232 
reduces K9 by one. Line 234 creates the A matrix by putting the matrix of exponents 
that was stored in B3 into BI. 

Lines 236-239 calculate the b vector. CC is a counter in the t’ subscripting 
system which controls the entry of b vector elements into B2(t’). Line 237 calculates 
the portion of the right hand side that is common to all terms. Line 238 multiplies that 
result by f(x) for objective function terms which Produces Py’ Ax’). Line 239 places 
InCp,- (x ")] into B2(t’). Lines 242-246 print rt on ‘) to the screen. Line 260 calls the 
ee linear equation subroutine which solves ALIn(y ‘y= =b. Lines 270-275 
print x to the screen and end the program. 


5. Subroutine To Stop Screen Printing, Figure 6.7 


600 INPUT"** Hit ENTER To Continue: "3Z9:RETURN 





Figure 6.7 Subroutine To Stop Screen Printing. 


Subroutine 600 is used to interrupt printing loops so that results are not 
scrolled off the screen before the operator can read them. 
6. Simultaneous Linear Equation Subroutine, Figure 6.8 
This subroutine solves simultaneous linear equations of the form Ax=b. K9 
is the dimension of the square A matrix and the x and b vectors. This subroutine 


documented in Appendix E, The Matrix Algebra Program. 


G. EXAMPLE PROBLEMS 
1. Example #1 
A design engineer wants to design a cylindrical oil storage tank with a storage 
capacity of 10007 cubic feet to put on an existing base. If the cost of construction 1s 


$1/foot? of tank surface, what are the optimal dimensions of the tank and how much 
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‘Simultaneous Linear Equation Subroutine: Ax=b 

‘Invert Matrix A 

FORK 7=K94+1TO2*K9: FORK8=1TOK9 
IFK7=K8+K9THENB1(K8 >K7 J=LELSEB1(K8 >K7 )=0 

NEXTK8:NEXTK7 

FORK7=1TOK9 

IFB1(K7>K7 )XSGN(B1(K7 >K7) )<1LE-8THENGOSUB 9910 

K2=1/81(K7 ,K7 ): FORK6=1TO2*K9:B1(K7,K6 )=B1IK7,K6 )¥K2:NEXTK6 
IFK7=K9THEN9865 

FORK8=K7+41TOK9: IFB1(K8>5K7 )=OTHEN9860 

K2=-B1(K8,K7 ) 
FORK6=K7TO2*K9:B1(K8 >K6 J=B1(K8,K6 )+(K2*B1(K7>K6 ) ) = NEXTK6 
NEXTK&8:NEXTK7 

FORK 7=K9TO2STEP-1 
FORK8=K7-1TOISTEP-1: IFB1(K8,K7 )=OTHEN9885 

K2=-B1(K8 >K7 ) 
FORK6=1TO2*K9:B1(K8 ,K6 J=B1( K8 >K6 )+(K2*B1(K7,K6 ) ) sNEXTK6 
NEXTK8: NEXTK7 

‘Mult A Inverse by b 
FORK7=1TOK9:B1(K7>1 )=0: FORK8=1TOK9:B1(K7,1 )=B1(K7,1)4+B1(K7 ,K8+K9 )*¥B2( K8 ) 
NEXTK8:NEXTK7: RETURN 

‘Error Routine 

IFERL>9700ANDERR=LITHENPRINT"!!!'ERROR: Matrix Is Not Invertable!!!";:END 
PRINT"Error Code";ERR3"In Line"3;ERL: END 

‘SWITCH ROWS 
FORK5=K7+41TOK9: IFB1(K5,K7 )*¥SGN(B1(K5 >K7 ) )<1LE-8THEN99G0 
FORK4=1TOK9*2: K3=B1(K7 >KG):B1(K7,KG )=B1L(KS KG ) 
B1(K5,K% J=K3: NEXTKG: RETURN 

NEXTK5:PRINT"Error: Matrix Not Invertable":END 





Figure 6.8 Simultaneous Linear Equation Subroutine. 


will it cost? The tank includes the cylindrical siding plus the top. The formulation 1s in 


Equations 6.26 and 6.27. 


Min $1(mr*) + $1(2mrh) (eqn 6.26) 
SUDJECL lo: 
mr-h = 1000% rh > 0 | (eqn 6.27) 


The objective function and constraint are posynomials, but the constraint is not in the 
< 1 form required by the program. Putting the constraint in S 1 form results in 
Equation 6.28. This is a zero degree of difficulty problem with three terms, two 


variables, one constraint, two terms in the objective function, and one term in the 


constraint. 


is 


1000r-h7! << 1 (eqn 6.28) 


The coefficients, c,’, 


t'=1,2,3 are m, 2m, and 1000 respectively. The input file and 
results of the program are at Figure 6.9. 


Input File: Results: 
o1 = lis. 55 = a bee 63 = 2/3 
1000 f(x) = $942. 48 


Optimal Values for r 





Figure 6.9 Input File and Results Of Example #1. 


The economic interpretation of 6, and 5 is that regardless of the price for 
steel, 1/3 of the cost will be for the top and 2/3 will be for the side. If the top and 
sides were constructed of different types of steel with different prices, these ratios 
would not change. 

2. Example #2 

This problem is the example stated in Figure 6.1 The input file is in Figure 
6.2 Dna, t’=1,2,3,4 are .5, .5, .5, and .75 respectively. The optimal value of the 
objective function is 40. The optimal values of erilas Ni7ay2) sabe 0, |, anda 


respectively. 
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Vil. MATRIX ALGEBRA PROGRAM 


A. GENERAL 

The matrix algebra program, MATALG, performs the following matrix algebra 
functions: matrix addition, multiplication, and inversion, scalar multiplication, 
calculation of determinants, integer exponentiation, and solutions to sets of 
simultaneous linear equations. MATALG is menu driven. The main menu enables the 
Operator to enter a new matrix, print the answer matrix, or call one of the functions 
listed above. Menus produced by each function prompt the operator for required 
input. Matrices may be entered from the M100 keyboard or from a RAM file. Output 
goes to the M100’s screen. Intermediate results may be displayed. Operations are 
performed in the conventional left to right order in which matrix operations are written 
out on paper,e.g. Ax Bx cl However, if the series of operations requires altering 
that order, the operator may store one matrix for future recall. Matrices may be 


entered in either the left or right hand position. 


B. TNPOT. 
The matrix input subroutine lets the operator select: 
©>> iiie position? of the matrix. If the new matrix is entered for the left side of the 
eeouen the program automatically places the old left side matrix on the right 
side. 


e Whether the matrix will be entered from the keyboard, the input file 
MATIN.DO, or from RAM storage. 


e Whether the matrix will be scalar multiplied or inverted. 

The matrix input routine must be accessed from the main menu to enter the first 
matrix. From then on, when a two matrix operation is selected from the main menu, 
the subroutine performing that operation automatically calls the input subroutine for 
the second matrix. 

If the input files MATIN.DO, is used, it must be created before the program 1s 
run. MATIN.DO may contain more than one matrix. Matrices must be preceded in 
MATIN.DO by their dimensions. An example of an input file is at Figure 7.1. 


See the general instructions on input files in Chapter 2. 


Left or right side of the operation. 


ce 


to in right contains two matrices: 


Figure 7.1 Sample Input File, MATIN.DO. 


When matrices are entered from the keyboard, the operator will be prompted for 


the matrix dimensions and for each matrix element. 


C. OUTPUT 

All output goes to the screen of the M100. Output may have up to three digits 
to the left and four digits to the right of the decimal point. If this configuration is not 
adequate, the operator may modify the format at the line numbers specified in Figure 


7.2 for the corresponding functions. 


LINE FUNCTION 
1082 Determinant. 
4075 Solution to Simultaneous Linear Equations 


6012 Other Matrix Output 
See the instructions for the PRINT USING command in Reference 1. 





Figure 7.2 Line Numbers Of Output Formats. 


D. DESCRIPTION OF VARIABLES 

e Al(3,K,K) holds the current matrices. Al(1 2 = Left side matrix/primary 
matrix/current. intermediate answer matrix. Al(2,,) = Right side/secondary 
matrix. Al(3,,) = Matrix being stored. 

e B1(K,K) holds the answer matrix as it is being calculated. 

e C(4) and Se, are the number of columns/rows in Al or Bl. &{3} and C3} 
correspond to oe ,). C2) and R(2) correspond to Al(Z,,), "Grand KG 
correspond to Al(3,,). C(4) and R(4) correspond to B1. 

e CC is the row counter in matrix output routine. 


¢ CD and RD are the column and row of element to be changed. 
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e CHI is the selection variable for the main menu. 
¢ DET(2) are the determinants of Al(1,,) and A1(2,,). 


e EF is the error Flag. 0.- No terminal error has been made. 1 - A terminal error 
has been made. A terminal error is one from which the program can not recover. 


e FF is the file Flag. 1 - MATIN.DO exists in RAM. 0 - MATIN.DO does not 
exist in RAM. 


ine. 45, 1, J2-andels are loop counters. 

¢ K is the largest dimension of largest matrix to be processed. 

e K4-K9 are counters in the simultaneous linear equation subroutine. 

e MEF 1s the multiplication/addition flag. 0 - Neither the multiplication nor the 
addition subroutines are running. | - The multiplication subroutine is running. 2 
- Addition subroutine is running. 


e MI is the matrix indicator. It shows which matrix in Al(1-3) is being operated 
upon. 


e MU is.an intermediate multiplier in the determinant and matrix inversion 
subroutines. 


e OF is the output flag. 1 - Send output to screen. 0 - Suppress output to screen. 


e SF is the simultaneous linear equation (SLE) flag: 1 - The SLE subroutine is 
running. 0 - The SLE subroutine is not running. 


e XP is the umber of times the matrix will be multiplied times itself in the integer 
exponentiation subroutine. | 


e Z9 is a general purpose variable. 


E. DESCRIPTION OF PROGRAM COMPONENTS 
A complete program listing 1s located at Appendix E. 


1. Initialization, Figure 7.3 


100 CLS:PRINT""': PRINT" *** MATRIX ALGEBRA PROGRAM ***";PRINT"" 
105 PRINT"IS INPUT MATRIX, 'MATIN.DO' IN RAM?": INPUT" OQO=NO, 1=YES"3FF 
107 IFFF=1THENOPEN"MATIN'FORINPUTAS1 


300 PRINT"**Enter The Single Largest Dimension of" 
305 INPUT"The Largest Matrix To Be Processed: "3K 
310 DIMA1(3>K,K),B1(K4+1,K¥2 ),R(4),0(74),DET( 2 ):MI=1:O0F=1:SF 





Figure 7.3 Initialization Section. 


Line 100 prints the program title. Line 105 permits the operator to specify 
whether file MATIN.DO will be used as a source of input and sets the file flag, FF, 
accordingly. Line 107 opens MATIN.DO for input if it 1s to be used. 


ue 


Lines 300-305 require the operator to specify the largest dimension, K, of the 
largest matrix to be processed. Line 310 dimensions matrices Al and BI and vectors 
R, C, and DET and initializes flags MI, OF, and SF. ° 

2. Main Menu, Figure 7.4 
This section permits the operator to select the next major operation to be 


conducted and calls the subroutine performing that operation. 


CLS: EF=0:PRINT"****MATRIX ALGEBRA PROGRAM MENU2xx«x«" 

PRINT" 1. Enter Starting Left Side Matrix" 

PRINT" 2. Matrix Inversion" 

PRINT" 3. Matrix Addition":PRINT" 4. Matrix Multiplication" 
PRINT" 5. Simultaneous Linear Equations" 

PRINT" 6. Print Current Answer Matrix":PRINT” 7. Other Options" 
INPUT" **Enter Number: "3CH 

IFCH=1LTHENMI=1:29=0:GOSUB7006 

IFCH=2THENMI=1:GOSUB2Z000 

IFCH=3THENGOSUB3000 

IFCH=4THENGOSUB5000 

IFCH=S5THENGOSUBGOO00 

IFCH=6THENGOSUB6000 

IFCH<>7THENGOTOSO1 

CLS: PRINT'"***MORE CHOICES:":PRINT” 1. Determinant” 

PRINT" 2. Matrix Integer Exponentiation" 

PRINT" 3. Store Current Matrix" 

PRINT" G. Retrieve Stored Matrix":PRINT" 5. Scalar Multiplication" 
PRINT" 6. Other Options": INPUT" **Enter Number: "3 CH 
IFCH=1THENMI=1:GOSUB1000 

IFCH=2THENMI=1 : GOSUB7600 

IFCH=3THENGOSUB8000 

IFCH=4THENMI=1 : GOSUB8200 

IFCH=STHENMI=1:GOSUB5100 

GNTOS501 





Figure 7.4 Main Menu. 


Line 501 initializes EF and prints the main menu header to the screen. Lines 
502-510 print the first screen of options and prompt the operator for a selection. Lines 
512-518 call the subroutine selected by the operator or branch to the second screen of 
options. Lines 520-532 print the second screen of options and prompt the operator for 
a selection. Lines 540-570 call the subroutine selected by the operator or return to the 
first screen of options. 

3. Pause Control Subroutine, Figure 7.5 
Lines 700-702 stop the program to permit the operator to view material on the 


screen and permit continuation by pressing the ENTER button. 
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"PAUSE CONTROL 


INPUT"** Hit ENTER To Continue"™3Z9:RETURN 





Figure 7.5 Pause Control Subroutine. 


4. Modifying the Secondary Matrix, Figure 7.6 


"INTERMEDIATE MODIFICATIONS 
PRINT"**Modify The 2nd Matrix?" 
INPUT" O=No, l=Invert, 2=Scalar Multiply: "329 


IFZ9=OTHENRE TURN 
MI=2: IFZ9=1THENGOSUB 2000ELSEGOSUB5100 
GOTO810 





Figure 7.6 Modifying The Secondary Input Matrix. 


This subroutine permits the operator to invert the second matrix of a two 
matrix operation or multiply that matrix by a scalar. Lines 810-812 print the options 
to the screen and prompt the operator for a selection. Line 815 causes a return 
without the matrix being modified if appropriate. Line 820 sets the matrix indicator, 
MI, to two and calls the matrix inversion or scalar multiplication subroutine. Line 825 
Starts the subroutine again, permitting the operator to select another option. 

5. Determinant Calculation, Figure 7.7 

Lines 1005-1008 test whether the matrix is square and print an error message 
if it is not. Line 1010 copies the matrix to be inverted into BI where the calculations 
will be conducted. Line 1020 initializes the value of the determinant as one and the 
row counter, I1. 

Line 1021 checks whether the diagonal element in the current row, R,, 1s zero. 
If so, then the row switching subroutine is called. If all the rows below R, have 0’s in 
column II then the determinant is zero. If a non-zero element can be found below R, 
in column I then that row is switched with R, and the determinant is multiplied by -1. 
Line 1022 tests for a terminal error from the row switching subroutine. Line 1023 


multiplies the determinant by diagonal element I1 and branches to the end of the 
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"CALC DETERMINANT 

IFR(MI )=C(MI JTHEN1010 

PRINT"ERROR: Number of rows/columns not equal:" 

PRINT" MATRIX IS NOT INVERTABLE !"*: GOSUB700:EF=1:RETURN 
FORI1=1TOR(MI ): FORIZ=1TOC(MI ):B1(1I1,I2 J=Al1(MI,I1,I12):NEXTI2:NEXTI1 
DET(MI )=1: FORI1=1TOR(MI ) 

IFB1(I1,I1)*¥SGN(B1(I1,I1) )<1LE-LOTHENGOSUB1 900ELSE1023 
IFEF=1THEN1008 

DET(MI J=DET(MI )*¥81(0I11,I1):IFI1L=R(MI JTHEN1O8O 

FORI3=1TOC(MI ):B1(1I1,13 )=B1(I1,13)/B81(I1,I1):NEXTI3 
FORI2=I1+1TOR(MI ): IFB1(12,I1 J=OTHEN1060 

FORI3=I1TOC(MI ):B1(1I2,13 )=B1(1I2,13)-(81(1I2,I1)*B1(1I1,13) ):NEXTI3 
NEXTI2:NEXTI1 

IFOF<>1THENRE TURN 

PRINT"**Det. Of Matrix ";MI3" Is: "5 
PRINTUSING"##FHE . RARE" 5 DET CMI ): GOSUB700 

RETURN . F 





Figure 7.7. Determinant Calculation. 


subroutine if the R, is the last row. Line 1025 divides all elements in R, by the 
diagonal element in R,. Lines 1030-1060 update the elements in R, and below in in 
column I1 and to the left. Line 1080 tests whether output is to be printed to the 
screen. If not, the subroutine ends. If so, lines 1081-1082 print the determinant. 


6. Row Switching Subroutine, Figure 7.8 


‘SWITCH ROWS 
FORJ=I1+1TOR(MI ): IFB1(J,I1)*SGN(B1(J,I1) )<1E-1OTHEN1940 
FORJ1=1TOC(MI )*¥2:TE=B1(I1,J1):B1(I1,J1 )=Bl(J,J1) 


B1(J,J1 )=TE :NEXTJ1:GOTO1950 
NEXTJ:EF=0:RETURN 
DET(MI J=-DET(MI ): RETURN 





Figure 7.8 Row Switching Subroutine. 


This subroutine is called when the determinant or inversion subroutines try to 
pivot on a row, Ryy, with a zero in the main diagonal element. The subroutine looks 
for the first row below Ry, that has a nonzero element in column 11.19 Line 1910 


searches the rows below Ry, for a row with a nonzero element in the appropriate 


10The same column as the zero on the main diagonal in Ry. 
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column.!! Lines 1920-1930 switch the elements of the rows using TE as an intermediate 
storage variable. If none of the rows below Ry, have a nonzero element in the 
appropriate column, then the matrix is not invertable and EF is set to one in line 1940. 
If a row switch was made, the determinant changes: sign in line 1950. 


7. Matrix Inversion Subroutine, Figure 7.9 


"MATRIX INVERSION 
OF =0:GOSUB1000: IFDET(MI )XSGN( DET (MI ) )>1E-1LOOREF=1THEN2017 
PRINT" ERROR: Determinant=0. MATRIX NOT INVERTABLE!":GOSUB700:EF=1 
IFEF=1THENRETURN 
FORI1=1TOR(MI ): FORI2=1TOC(MI ):B1(I1,I2)=AL(MI ,I1,12):NEXTI2:NEXTI1 
FORI1=C(MI J+1TO2*C(MI J: FORIZ=1TOR(MI ) 
IFI1=I2+R(MI JTHENB1(12,I1 J=1LELSEB1(I2,I1)=0 
NEXTI2:NEXTI1L 
FORI1=1TOC(MI ) 
IFB1(I1,I1)*SGN(B1(1I1,I1))<1E-1LOTHENGOSUB1900ELSE2055 
IFEF=1THEN2015 
MU=1/B1(1I1,I1):FORIZ=1TO2xC(MI ):B1(I1,13 J=B1(1I1,13 JXMU:NEXTI3 
IFI1=C(MI JTHENZ080 
FORI2=I1+17TOR(MI J): IFB1(12,11)=OTHEN2075 
MU=-B1(1I2,I1) 
FORI3=I1TO2*C(MI ):B1(1I2,13 J=B1( 12,13 )+(MUXB1(11,13)):NEXTI3 
NEAXATI2:NEXTI1 
FORI1=C(MI JTO2STEP-1 
FORI2=11-17TO1STEP-1: IFB1( 12,11 )=OTHEN2130 

=-B1l(I2,I1) 
FORIZ=1TO2*C(MI ):81(12,13 J=B1( 12,13 )+(MUXB1(1I1,I3) ):NEXTI3Z 
NEXTI2:NEATI1 
FORII=1TOC(MI J): FORIZ=1TOR(MI ) 
A1(MI,I2,I11)=81(12,I1+C(MI) JsNEXTI2Z:NEXTI1 
OF=1:RETURN 





Figure 7.9 Matrix Inversion Subroutine. 


The matrix inversion subroutine places the matrix to be inverted, ft, and an 
identity matrix, I, into BI. Each row of BI holds a row of wt and the corresponding 
row from |. Elementary row operations are conducted on pt in BI to change that 
portion of BI to an identity matrix. The same elementary row operations are 
conducted on the portion of BI that started as an identity matrix. When the portion 
of BI which started as pt is changed to I, then the portion of BI which started as | 


becomes wl 


l1The decision rule actually looks for an element outside the range 0 + 19°10. 
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Line 2010 stops intermediate results from being printed to the screen by 
setting OF to zero. Line 2010 alsoealismhe determinant calculation subroutine and 
tests whether the determinant is equal to zero. If the determinant equals zero, then 
lines 2015-2017 print an error message, set the error flag, and terminate the inversion 
subroutine. Line 2020 copies pt into Bl. Lines 2030-2035 place an identity matrix with 
the same dimensions as pt into B1 with pp. 

Lines 2040-2075 Line 2045 checks whether the MrT] is zero. If so, then the 
row switching subroutine is called. If the row switching subroutine can not find a row 
for which element I] does not equal zero, then the matrix is not invertable and line 
2046 branches to the error message. Line 2055 calculates the constant, MuU,!? and 
multiplies row Il by MU. Since lines 2060-2075 do not apply to the last row of ph, line 
2057 branches around them if I] points to the last row. 

Lines 2060-2075 perform the elementary row operations to change to zero the 
elements of column I] that are in rows below row il. Lines 2080-2130 perform the 
elementary row operations which change to zero the elements of ft above the main 
diagonal. Lines 2140-2145 copy the inverted matrix from BI back to the appropriate 
section of Al. Line 2190 turns the output back on by resetting OF and terminates the 
subroutine with a return. 

8. Matrix Addition, Figure 7.10 


3000 ‘MATRIX ADDITION 
3010 MF=2:G0OSUB 7000: GOSUB800: IFEF=1THENRETURN 


3015 FORI1=1TOR(1):FORI2=1TOC(1):A1(1,I1,12 J=A1(1,11,12)+A1(2,11,12) 
3020 NEXTI2:NEXTI1:GOSUB6000 : MF=0: RETURN 





Figure 7.10 Matrix Addition Subroutine. 


Line 3010: 


¢ Sets MI=2 indicating to the input subroutine that it is being called from the 
matrix addition subroutine. 


e Calls subroutine 7000 to enter the second matrix. 


a. 12\fultiplying row Il by MU makes main diagonal element Mr equal to one, 
1.e, its identity matrix value. 
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e Calls subroutine 800 to permit the operator to invert the second matrix or 
multiply it by a scalar. 


e Evaluates whether a terminal error was made in either subroutine 7000 or 800 
and, if so, terminates the matrix addition subroutine. 


Lines 3015-3020 add the elements of Al(1,,) and A2(2,,). Line 3020 also calls 
subroutine 6000, printing the answer, resets MF=0, and terminates the matrix addition 
subroutine. 


9. Simultaneous Linear Equations, Figure 7.11 


‘SIMULTANEOUS LINEAR EQUATIONS 

CLS:PRINT"*¥*¥Solves Ax=b. Choices:":PRINT" 1. Enter b Vector" 
PRINT" 2. Change An Element In Matrix A" 

PRINT" 3. Solve Current Ax=b" 


e 


PRINT" 4. Return": INPUT" * Select A Number: '"3CC 
IFCC=1G0T04040 

IFCC=2GO0T04050 

IFCC=3G0T04060 

IFCC=GTHEN RETURN 


GOTO 4000 

MI=2:R(2)=C(1):C( 2 J=1;GOSUB7040: GOTO04000 

INPUT"**Row, Column Of Matrix A To Be Changed: "3RD,CD 

PRINT" - Enter Row"5$RD3", Column"53CD3":"3:INPUTA1(1,RD,CD):GOTO4000 
MI=1:SF=1:0F =0:GOSUB8000: GOSUB2000 

IFEF=OTHENGO70 

PRINT"*Solution Not Uniquely Determinable" :GOSUB700: RETURN 
GOSUB5000 :CC=0: FORI1=1TOR( 2 ):CC=CC+1:PRINT"x("3I135") = "3 
PRINTUSING' FHF. HEHEHE" 3B1(1I1,1):IFCC>6THENGOSUB700:CC=0 
NEXTI1:GOSUB700:SF=0: GOSUB8200:GOT04000 





Figure 7.11 Simultaneous Linear Equation Solving Subroutine. 


This subroutine solves sets of linear equations of the form Ax=b where A 1s 
m by m matrix of rank m and x and b are vectors of length m. Matrix A must be 
entered as the primary matrix before this subroutine is called. The subroutine prompts 
the operator to enter b. The subroutine also permits the operator to change individual 
elements of the A matrix. 

Lines 4010-4014 print a header and a menu of options and prompt the 
operator to select an option. Lines 4020-4035 transfer control to execute the option 
selected. Line 4040 sets MI =2, indicating that b will be stored in A1(2,,), dimensions 
the b vector, and calls subroutine 7040 to input b. Line 4050 prompts the operator for 
the row and column of the element in A to be changed. Line 4052 prompts the 


operator to enter the new value for that element. 
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Lines 4060-4080 solve the system of equations. Line 4060 sets MI, SF, and 
’ OF and calls subroutines which store, then invert, the A matrix. Lines 4064-4065 print 
an error message if the A matrix is not invertable. Lines 4070-480: 
¢ Multiply Av! by b producing the solution vector, x. 
e Print the solution vector. 
e Reset SF=0. 
e Retrieve the stored A matrix. 


10. Matrix Multiplication Subroutine, Figure 7.12 


‘MATRIX MULT 

MF=1:IF SF=1THENS5020 ‘ 

MI =2:GOSUB 7000 : GOSUB800 : [FEF=1THENRETURN 

R(4)=RE1 ):C0064)=C(O 2): FORTL=1TOR(4 ): FORT 2=1T0C(4):B1l(I1,12)=0 


FORI3=1TOC(1):B1(I1,12 )=A1(1,11,13 )*A1(2,13,12)+B1(I1,12) 
NEXTI3:NEXTI2:NEXTI1:MF=0 

IF SF=OTHENGOSUS37500:GOSUB6000 

RETURN 





Figure 7.12 Matrix Multiplication Subroutine. 


Line 5010 sets MF!? and branches to avoid the input subroutines in line 
5015 if SF equals one.!4 Line 5015 calls the subroutines to enter and modify the second 
matrix. Lines 5020-5024 set the dimensions Bl and perform the multiplications and 
additions required to place the answer matrix in BI. If SF equals zero,!> then line 
5050 calls subroutines which copy the answer from BI to Al(1,,) and print the answer 
matrix. 

11. Scalar Multiplication Subroutine, Figure 7.13 
This subroutine multiplies matrix AI(MI,,) by a scalar, SM. Lines 


5110-5115 prompt the operator to enter the scalar and conduct the multiplication. 


I3\{F =] indicates that matrix multiplication is being performed. 


l4That is, if the matrix multiplication subroutine is called from the simultaneous 
linear equation subroutine. 


1>That is, this subroutine is not being called from the simultaneous linear 
equation subroutine. 
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5100 'SCALER MULT 


5110 INPUT"Enter Scalar Multiplier:";SM:FORI1=1TOR(MI ): FORI2Z=1TOC(MI ) 
5115 Al(MI,I1,I2)=AL(MI,I1,12 JXSM:NEXTI2:NEXTI1: RETURN 





Figure 7.13 Scalar Multiplication Subroutine. 


12. Subroutine To Print Al(1,,), Figure 7.14 


‘PRINT OUTPUT MATRIX 
PRINT" ** Current Answer Matrix:':CC=0:FORIL=1LTOR(1):CC=CC+l1 


FORI2=1TOC(1):PRINTUSING' H#H# . HHH" 5A1(01,11,12)3 sNEXTI2:PRINT"" 
IFCC=3THENGOSUB 700 :CC=0 
NEXTI1:GOSUB700 : RETURN 





Figure 7.14 Subroutine To Print The Primary Matrix. 


Lines 6010-6070 print a header and then print Al(1,,). The matrix is printed 
up to three rows at a time. 
13. Matrix Input Subroutine, Lines 7000-7050 
a. Input Matrix Configuration, Figure 7.15 
The section in Figure 7.15 prompts the operator to specify: 


e Whether the matrix to be entered will go on the left or right hand side of the 
Operation. 


e Whether the matrix will be entered from the keyboard, MATIN.DO, or retrieved 
from RAM storage. The dimensions of the incoming matrix are entered from the 
appropriate source. 

Line 7001 prompts the operator to specify whether the incoming matrix will 
be on the left or right side of the operation. If the incoming matrix is to be on the left 
side, lines 7003-7004 move the matrix in Al(1,,) to Al(2,,). Lines 7006-7008 print an 
appropriate header. Lines 7009-7011 print the source options for the incoming matrix. 
The option to enter a matrix from MATIN.DO will be printed only if MATIN.DO has 
been created (FF=1) and the end of MATIN.DO has not been reached (EOF(1)=0). 
Lines 7012-7013 transfer control to enter the matrix from the appropriate source. 


Lines 7014-7018 enter the dimensions of the new matrix from the appropriate source. 
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‘MATRIX INPUT 

CLS:PRINT"":PRINT"Will This Matrix Be On:":INPUT” O=Left, 1=Right"3Z9 
IFZ9=1THEN7006 

R(2)=R(1):C( 2 )=C( 1): FORI1=1TOR( 1): FORI2=1TOC(1):A1(2,11,I12)=A1(1,1I1,1I2) 
NEXTI2:NEXTI1:MI=1 2 

CLS: IFMI=2THEN7008 

PRINT'"**Choices For Left Hand Matrix':GOTO7009 

PRINT"xxChoices For Right Hand Matrix:" 

PRINT" 1. Enter Matrix From Keyboard" 

PRINT" 2. Retrieve Stored Matrix": IFFF<>1THEN7012 

IFEOF(1)=OTHENPRINT" 3. Enter Matrix From MATIN.DO" 

INPUT" xxEnter A Number: "3Z9:IFZ9=1THEN7015 

IFZ9=3THEN7018 

R(MI J=R(3): COMI )=C(3 ):GOTO7020 

PRINT" **Enter The Rows, Columns" 

INPUT"In The Next Matrix: "“3;R(MI),C(MI ):GOTO7020 

INPUT#1,R(MI),C(MI ) : ; 





Figure 7.15 Input Matrix Configuration. 


b. Detection Of Dimensioning Errors, Figure 7.16 


IF MF<>1THEN7030 

IFR( 2 )=C(1JTHEN7030 

PRINT'**ERROR: Columns in LEFT MATRIX ="3C(1) 
PRINT" Rows In Right Matrix ="5R(2) 


PRINT"These Must Be Equal For Matrix Mult! !":GOSUB700:EF=1:GOTO7006 
IF MF<>2THEN7035 

IF(R(C1)=RO2 ANDC(1)=C(2)) THEN 7035 

PRINT'*x*XERROR:Dimensions For Both Input" 

PRINT'Matrices Must Be Equal!t":GOSUB700:EF=1:GOTO7006 





Figure 7.16 Detection Of Dimensioning Errors. 


If the matrix being entered is the second matrix in a matrix multiplication 
operation, then lines 7020-7026 check whether the number of columns in the left matrix 
is equal to the number of rows in the right matrix. If not; then an error message is 
printed and control is transfered to the beginning of the matrix input subroutine. If 
the matrix being entered is the second matrix in a matrix addition operation, then lines 
7030-7034 check whether the dimensions of the left and right matrices are the same. If 
not, then an error message is printed and control is transfered to the beginning of the 


matrix input subroutine. 
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c. Matrix Input Section, Figure 7.17 


7035 IfZ9=2THENGOSUB8200: RETURN 

7036 IFZ9=3THEN7050 

7037 PRINT" *¥Fill Matrix Row By Row:":PRINT"" 

7040 FORIL=1TOR(MI J): FORIZ=1TOC(MI ) 

7042 PRINT"-Enter Row';I13;"And Column"3I23":"3 

7044 INPUTAIL(MI,I1,12):NEXTI2:PRINT'":NEXTI1: RETURN 

7050 FORII=1TOR(MI ): FORIZ=1TOC(MI ): INPUT#1,A1(MI,I1,12):NEXTI2:NEXTI1: RETURN 








Figure 7.1/7 WNlatrix Input Section. 


If the incoming matrix is to be retrieved from RAM storage, line 7035 calls 
the appropriate subroutine. Line 7037-7044 enter the incoming matrix from the 
keyboard. If the incoming matrix is to be entered from MATIN.DO then line 7036 
transfers control to 7050 which performs the entry. 

14. Copy B1 Into A1(1,,), Figure 7.18 


7500 ‘COPY Bl INTO Al(l,;,) | 
7510 R(1)=R(04):C01)=C(4): FORI1=1TOR( 1): FORIZ2=1T0C(1) 
7512 Al(1,I1,12)=B1l(I1,12):NEXTIZ2:NEXTI1: RETURN 


Figure 7.18 Subroutine to copy B1 into Al(1,,). 


Lines 7510-7512 dimension al(1,,) and copy B1 into A1(1,,). 


15. Matrix Integer Exponentiation, Figure 7.19 


7600 "MATRIX INTEGER EXPONENTIATION 
7610 CLS: PRINT": INPUT"**Enter Integer Exponent > 2: '" 5XP 


7620 R(2)=R(1):C0 2 )=C( 1): FORIL=1TOR( 1): FORI2Z=1TO0C( 2) 
7622 41(2,11,12)=A1l(1,I1,12):NEXTI2:NEXTI1 
7630 SF=1: FOREX=2TOXP : GOSUB5020 : GOSUB7500 : NEXTEX: GOSUB6000 : SF=0: RETURN 





Figure 7.19 Matrix Integer Exponentiation Subroutine. 
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This subroutine raises the primary matrix to an integer power greater than or 
equal to two. Line 7610 prompts the operator to enter an exponent, XP. Lines 
7620-7622 dimension Al(2,,) and copy Al(1,,) into A2(2,,). Line 7630: 

© Sets SF=1. This suppresses printing of intermediate results. 
e Performs XP-1 matrix multiplications. 
e Prints the final result. 


16. Storage and Retrieval Subroutines, Figure 7.20 


"STORE Al(1,,) 
R(3)=R(1):C03 J=C(1): FORIIT=1TOR(3 ): FOREZ=1T0C(3 ) 
A1(3,I1,I12)=A1(1,I1,I2):NEXTIZ:NEXTI1: RETURN 


"RETRIEVE THE STORED MATRIX 
RUMI J=R(3):COMI J=C(3 ): FORIL=1TOR(MI ): FORIZ=1TOC(MI ) 
A1(MI ,I1,I2)J=A103,I1,12):NEXTIZ2:NEXTI1: RETURN 





Figure 7.20 Storage and Retrieval Subroutines. 


Lines 8010-8012 dimension A1(3,,) and store Al(1,,) in A1(3,,). Lines 
8210-8212 dimension Al(1,,) and retrieve Al(1,,) from A1(3,,). 


F. SIMULTANEOUS LINEAR EQUATION SUBROUTINE, FIGURE 7.21 

Many programs require a simultaneous linear equation solver. Often these 
programs compute the A matrix as part of the program and use the results in 
subsequent calculations. The following subroutine may be inserted in other programs 
without requiring the loading of the entire matrix algebra prograin. 

This subroutine follows the same algorithm as the simultaneous linear equation 
subroutine in the Matrix Algebra Program. K9 is the dimension of the A matrix. Bl 
holds the A matrix; B2 holds the b vector. The x vector is returned in the first column 
of Bl. If the A matrix is to be used later, it must be stored somewhere other than B1 
since the A matrix in B1 is changed to an identity matrix by this subroutine. Instead 
of testing for a zero determinant, the subroutine uses the error identification subroutine 
9900 to determine if the A matrix is not invertable. Variables K2-K9 are used to avoid 


conflict with other counters. 
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"Simultaneous Linear Equation Subroutine: Ax=b 

DIM B1(K941,K9¥2),B2(K9);'Bl = A matrix; B2 = b vector 
"Input from SLEIN.DO; Set MAXFILES in main program 
OPEN" SLEIN"FORINPUTAS9 

FORK8=1LTOK9: FORK7=1TOK9: INPUT#9,B1(K8,K7):NEXTK7:NEXTK8 
FORK8=1LTOK9: INPUT#9,B2(K8 ):NEXTK8 

‘Invert Matrix A 

FORK 7=K94+1TO2*K9: FORK8=1LTOK9 

IFK7=K8+K 9THENB1( K8,K7 )=LELSEB1( K8,K7 )=0 

NEXTK8: NEXTK7 

FORK 7=1LTOK9 

IFB1(K7,K7 )¥SGN(B1(K7,K7) )<LE-8STHENGOSUB9910 
FORK6=1LTOZ2*K9:B1(K7 ,K6 J=BL(K7 ,K6 )/BIL(K75K7) sNEXTK6 
IFK7=K9THEN9865 

FORK8=K 741 TOK9: IFB1(K8,K7 )=OTHEN9860 

FORK6=K 7TO2*K9:81(K8,K6 J=B1(K8,K6 J-(B1(K8,K7 )¥B1(K7>,K6) ) sNEXTK6 
NEXTK8 : NEXTK7 , 
FORK 7=K9TO2STEP=-1 
FORK8=K7-1TOISTEP=-1: IFB1(K8,K7 )=OTHEN9885 
FORK6=1TO2*K9:B1(K8 ,K6 J=B1(K8,K6 J-(B1(K8 yK7 )¥B1(K7,K6) ) sSNEXTK6 
NEXTK8: NEXTK7 

"Mult A Inverse by b 

PRINT"** Sim Lin Eq Solution: X(i) =" 

FORK 7=1LTOK9:B1(K751 J=1: FORK8=K941TOZ2*K9: FORK6=1LTOK9 
B1(K7,1)=B1(K7>,1)4+B1(K7,K8 )*82(0K6 ): NEXTK6:NEXTK8 
PRINTUSING "#8384. ###"5B1(K7,1):NEXTK7: PRINT": RETURN 
‘Error Routine 

IFERL>97Q0QANDERR=LITHENPRINT"!! "ERROR: Matrix Is Not Invertable!!!";END 
PRINT"“Error Code";ERR3"In Line"3;ERL:END 

"SHITCH ROWS 
FORK5=K7+1LTOK9: IFBIL(KS,K7 )¥SGN(B1(K5,K7) )<LE-8THEN9940 
FORKG=1TOK9¥2 :K3=B1(K7 KS ):B1(K7 KS )=B1(K55,K4 ) 

B1(K5,K4 )=K3:NEXTKG: RETURN 

NEXTK5:PRINT"Error: Matrix Not Invertable":END 





Figure 7.21 Simultaneous Linear Equation Subroutine. 
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VUT. NUMERICAL DOUBLE INTEGRATION PROGRAM 


A. GENERAL 

This program numerically integrates functions of one or two variables using 
Simpson’s Rule with a Romberg extrapolation to improve accuracy. The Romberg 
extrapolation is described in [Ref. 12:pp.250-276]. Using the Romberg extrapolations 
allows the operator to specify an acceptable error. The program conducts 
extrapolations until the error of the numerical estimate is below that specified 
tolerance. 

The operator may interactively change the function being integrated, the limits of 
integration, or the Romberg tolerance. This program is also written as a subroutine. 
However, in the subroutine the function being integrated, the limits of integration, and 


the tolerance may not be interactively changed. 


B. INPUT 

All input is entered from the keyboard of the M100. When the program begins, 
a menu appears which allows the operator to select whether the function to be 
integrated, the limits of integration, or the Romberg tolerance will be changed. 

When the operator selects an input to be changed, the program calls the edit 
function for the applicable lines in the program. The edit function terminates the 
running of the program. The operator should change only the right hand side of the 
input equations. After changing the lines required, the operator will hit the F8 button 
on the M100. This puts the M100 back in the BASIC mode. The operator must then 
enter RUN or hit the F4 button to run the program again. If the operator wants to 
change another input, he should select another input from the menu and repeat the 
process. 

1, Changing The Function, f(x,y), To Be Integrated 

Enter zero at the main menu. When lines 1285-1288 appear, change the right 
hand side of the equation on line 1286. If the equation is too long for one line, then: 
¢ Calculate a partial function value on line 1286, assigning it to F. 


e Add a line 1287 ASIAN: the final function value to_F and including the partial 
function value from line [286 in the right hand side of the equation on line 1287. 


a2 


2x3 
For example, if f(x,y) = (x2 +y°+7) * eX “Y', then the function might be broken 


down as indicated in Figure 8.1. 


1286 F=EXP(X* 2xy%3 ) 


1287 FH(X*2+Y%5+7 )xF 





Figure 8.1 Example Of Function To Be Integrated. 


Up to two independent variables, X and Y, may be used in the equation. The operator 
must ensure that f(x,y) is formulated with X as a variable for which constant limits of 
integration can be specified. After f(x,y) is entered, depress F8, then F4 to return to 
the main menu. 
2. Changing The Limits Of Integration 
Enter one at the main menu. When lines 1291-1298 appear, change the right 
hand side of the equations on lines 1293-1298 as desired. The upper and lower limits 
of integration for X, XUPPER and XLOWER, must be constants. The upper and 
lower limits of integration for Y, YUPPER and YLOWER, may be constants or given 
in terms of X. Do not alter the return statement at line 1295. After the limits of 
integration are entered, depress F8, then F4 to return to the main menu. 
3. Changing the Romberg Tolerance 
The operator should enter two from the main menu and enter the new 
tolerance when prompted. 
4. Using The Program For Single Integration 
Although the program is written for double integration, single integration may 
be calculated using the following steps. 
e Set the function to be integrated equal to one, te. line 1286 will be F=1. 
¢ In lines 1296-1297 set YUPPER=f(X) and YLOWER=0. 
e In lines 1293-1294 set XUPPER and XLOWER as the constant limits of X 
between which the function YUPPER= f(X) 1s to be integrated. 


Porexamples tor i, x? dx, the corresponding limits of integration in lines 1293-1297 


would be as indicated in Figure 8.2. 


1293 XUPPER=2 
1294 XLOWER=0 


1295 RETURN 
1296 YUPPER=X*% 2 
1297 YLOWER=0 





Figure 8.2 Example of Limits Of Integration. 


C. OUTPUT 

The estimated value of the integral is printed to the screen with its tolerance 
error. If the program generated a tolerance error that was less than the tolerance 
specified 1n the input, then that tolerance is printed. If the program could not generate 
an estimate within the specified tolerance, then a message to that effect is printed to 


the screen. 


D. EXPLANATION OF VARIABLES 

¢ A2(6,6) is the matrix holding Romberg extrapolation values. 

e DX and DY are the widths of intervals (XU-XL)/N and (YU-YL)/N respectively. 

e F is the value of f(x,y) to be integrated at a particular point. 

e Ji through J9 are loop counters. 

e N is the number of intervals into which the distances XU-XL and YU-YL are 
divided. 

e SS is the Simpson’s Rule sum specified in equation 8.1. In equation 8.1 
fee f(x, ylx= X), YLSy;s YU, 1=1,2,3,....n+1. mis the number of intervals into 


which the distance YU-YL has been divided. n must be a positive, even integer. 
Simpsons’s Rule Sum = f, +4f,+2f,+4f,+ 2f5+,...,+4f +f 4 1 (eqn 8.1) 


e TL is the user specified tolerance. 

¢ XLOWER or XL is the lower limit of integration for X. 
e XUPPER or XU 1s the upper limit of integration for X. 
e YLOWER or YL 1s the lower limit of integration for Y. 
e¢ YUPPER or YU is the upper limit of integration for Y. 


e Z9 1s a Selection variable. 
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FE. EXPLANATION OF PROGRAM COMPONENTS 
A complete program listing is located at Appendix F. 


1. Initialization, Figure 8.3 


1200 'Numberical Double Integration:Steven H. Cary:24¢ Aug 86 


1201 DIMA2(6,6):TL=.001 





Figure 8.3 Initialization Section. 


Line 1201 dimensions the matrix holding the Romberg extrapolations and sets 
the default tolerance to .001. 
2. Option Selection, Figure 8.4 


CLS:PRINT’": PRINT" ** Double Integration %*%*" 
PRINT" Romberg Algorithm" 

PRINT"O=Edit Function To Be Integrated." 

PRINT" 1=Edit Limits Of Integraticn.” 


PRINT" 2=Edit Tolerance; Current Tol.="3TL 

PRINT" 3=Calculate Integral ":INPUT"Enter 0, 1, 2,5) or 3:"32Z9 
IFZ9=OTHENEDIT1285~-1288 

IFZ9=1THENEDIT1291-1298 
IFZ9=2THENPRINT"": INPUT" Tolerance="3TL:GOTO1205 





Figure 8.4 Option Selection Section. 


Lines 1205-1218 print a menu which permits the operator to change f(x,y), the 
limits of integration, or the tolerance, or to calculate the integral. If f(x,y) or the limits 
of integration are to be changed, then lines 1216 or 1217 activate the editor for the 
appropriate program lines. The editor terminates program execution, thereby requiring 
that the program be executed after editing. If the tolerance is to be changed, then line 
1218 prompts the operator to update TL and redisplays the menu. 

3. Integration Calculation, Figure 8.5 

Line 1220 clears the screen during the calculation and prints an admonition to 

be patient while the calculation occures. Line 1230 sets the initial number of intervals 


to two, calls subroutine 1293, which calculates the interval width, DX. Line 1240 


a5 


CLS: PRINT": PRINT" !'!Be Patient!!":PRINT"" 
N=2:GOSUB1293 : DX=(XU-XL )72 

FORJ9=1T06 : DX=DX/2 s N=N¥2 ; 

X=XU : GOSUB1296 : GOSUB1280 :A2(J9»1 J=SS*DY 


X=XL : GOSUB1296 : GOSUB1280 : A2(J9>1 J=A2(J9,1)4+SS*DY 
FORJ8=2TON: X=X+DX: GOSUB 1296 : GOSUB1280 
A20J59,1I=A20 59,1 J+ 2¥SSHDY : NEXT JUS 
A2(J9,1)=A2( 59,1 *XDX/3 





Figure 8.5 Integration Calculation. 


starts a loop in which the Simpson’s Rule intervals are halved at each iteration. That. 
is, in the first iteration YU-YL and XU-XL are divided into four intervals, in the 
second iteration they are divided into eight intervals, and so on for six iterations. Line 
1240 cuts the interval for X, DX, in half and doubles the number of intervals, N. 

Lines 1242-1245 call the subroutines which compute the Simpson Rule sums, 
SS, at the upper and lower bounds of X. These sums are multiplied by their respective 
interval widths, DY, and added together into A2(J9,1). Lines 1250-1251 calculate the 
same summation for values of X between XL and XU at intervals DX and and add the 
sums to A2(J9,1). Line 1252 multiplies A2(J9,1) by DX/3 to complete the Simpson’s 
Rule approximation of {yj f(x,y) dydx. 

4. Romberg Extrapolation, Figure 8.6 


1255 IFJ9=1THENNEXTJ9 


1260 FORJ8=1TOJ9-1 
1261 A2(J9,J84+1 )=A2( 59,58 J4( (A205 9 »J8 J-A2( 59-1558 ) 714% J8-1) )sNEXTIS 





Figure 8.6 Romberg Extrapolation. 


Because the Romberg extrapolations require two numerical approximations, 
line 1255 skips the extrapolation section after the first iteration, i.e. when J9=1. Lines 
1260-1261 conduct the Romberg extrapolation as described in the section on Romberg 
extrapolation in [Ref. 12:pp. 250-276]. 
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5. Termination and Output, Figure 8.7 


T1=A2059,J59)-A2(59,J9-1 ): IFSGN( T1 )*¥T1-TLOOTHENNEXTJ9ELSE1264 
PRINT’ Tolerance of";TL3"not met after five extrapolations" 
IN=A2(J9,J9) 

PRINT" Integral ="3:PRINTUSINGHH##H8H . HHHHH" SIN 

PRINT” Actual Tolerance=";:PRINTUSING''##. ##4#" 3 T1L*XSGN(TL) 


SOUND1567,10:SOUND1244 ,10: SOUND1046 ,10: SOUND 783 ,20 
SOUND1046 > 10: SOUND 783 »40 

INPUT“Hit Enter To Continue: "3Z9:GO0TO1205 
FORJ7=1T06: FORJ6=1TOJ7: PRINTUSING"##. ###" 5A2(J7,J6)3 
NEXTJ6: PRINT" :NEXTJ7: INPUTZ9: RETURN 





Figure 8.7. Program Termination and Output. 


Line 1262 finds the difference, Tl, between the last two Romberg 
extrapolations and compares that difference to the user specified tolerance, TL. If the 
difference is greater than the tolerance, then the extrapolation is not accurate enough, 
another numerical integration is conducted with DX and DY halved, and another 
extrapolation is made. Otherwise, the integral is within tolerance and the program 
branches to line 1264 to begin the output sequence. If the integral is not within 
tolerance after six iterations, iterations,!® the program terminates. Line 1263 prints a 
message to the screen indicating that the tolerance has not been met. Line 1264 
assigns the final value of the integral to IN. Lines 1265 and 1266 print the value of the 
integral and the actual tolerance,!’ to the screen. Lines 1267 and 1268 play a short 
tune to cue the operator that the calculation has finished. Line 1269 holds the results 
on the screen until the operator hits ENTER, cycling the program back to the main 
menu at line 1205. 

6. Diagnostic Subroutine, Figure 8.8 

Lines 1275-1276 print the A2 matrix. This subroutine can be called in the 
middle of a calculation to check how far the calculation has progressed. To call the 
subroutine in the middle of a calculation: 

e Hit SHIFT and BREAK together to stop the calculation. 
Enter GGounl2 7/5. 


loThat is, after distances XU-XL and YU-YL have been broken into 128 
intervals. 


17 Actual tolerance may be less than the user specified tolerance. 


on 


1275 FORJ7=1T06: FORJ6=1TOJ7: PRINTUSING'## . ###" 5A2( 57 »J6 ) 


1276 NEXTJ6:PRINT'":NEXTJ7: INPUTZ9: RETURN 





Figure 8.8 Diagnostic Subroutine, Prints Matrix A2?.. 


e After viewing the matrix, hit ENTER to continue the program. 


7. Simpson’s Rule Summation, Figure 8.9 


1280 REM Simpson's Rule Sum 

1281 Y=YU:GOSUB1285:SS=F :Y=YL:GOSUB1285:SS= SS+F 

1282 FORJ5=2TO(N/2 ): Y=Y+DY : GOSUB1285 : SS=SS+4*F : Y=Y+DY : GOSUB1285 
1283 SS=SS+2*F :NEXTJ5: Y=Y+DY : GOSUB1285: SS=SS+4*F : RETURN 





Figure 8.9 Simpson’s Rule Summation Subroutine. 


Lines 1281-1283 calculate SS= )) Ff, + 4f + 2f, + 4fg + 2f— + .....4f, + 
eri where f= f(x,ylx=X), Nils ys UE i= 1,2,3,...nt+1. m is the number of 
intervals into which the distance YU-YL has been divided. 


8. F(x,y) to be integrated, Figure 8.10 


‘f(x, ,y) to be integrated: 
F=1 


'X & Y=independent variables. Hit F8, then F4 When Done. 
RETURN 





Figure 8.10 Subroutine To Calculate f(x,y). 


The function to be integrated, f{x,y) 1s at line 1286.18 Lines 1285 and 1288 are 


comments printed to the screen during editing to assist the operator. 


I8An additional line, 1287, may be added if the function is too long for one line. 
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9. Limits Of Integration, Figure 8.11 


‘Limits of Integration: 

"XLOWNER/XUPPER are constants. 

‘YUPPER & YLONER may be constants or given in terms of X. 
XUPPER=1.5707963 

*XLOWNER=0 


RETURN 

YUPPER=SIN(X ) 

YLOWNER=0 

"Hit F8, Then F4 When Done 
DY=(YU~-YL )/(N+1): RETURN 





Figure 8.11 Limits Of Integration Subroutines. 


Upper and lower limits of integration for X are entered at lines 1293 and 1294 
respectively and must be constants. Limits of integration for Y are entered at lines 
1296-1297 and may be either constants or functions of X. Line 1299 updates DY. 


Lines 1290-1292 and 1298 are comments to assist the operator during editing. 


F. INTEGRATION SUBROUTINE 

The numerical integration program described above is adapted in Figure 8.12 for 
use aS a subroutine. In the subroutine neither f(x,y), the limits of integration, nor the 
tolerance can not be edited during program execution. All comment lines to facilitate 
editing have been removed. The subroutine returns IN as the numerical approximation 
of the integral but does not print IN. The operator must dimension A2(6,6) with the 


other arrays in the main program and delete line 1201 in the subroutine. 
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1200 
1201 
1202 
1220 
1230 
1240 
1242 
1245 
1250 
1251 
1252 
1255 
1260 
1262 
1263 
1264 
1266 
1275 
1276 
1280 
1281 
1282 
1283 
1286 
1289 
1293 
1294 
1295 
1296 
1297 
1279 


‘Numberical Integration Subroutine:Steven H. Cary:24 Apr 86 
DIMA2(6,6) 

TL=.001 

CLS:PRINT"''; PRINT" ''Calculating An Integral!!":PRINT"’” 
N=2:GOSUB1293 :DX=( XU-XL )72 

FORJ 9=1T06 : DX=DX/2 : N=N¥2 

X=XU ;: GOSUB1296 : GOSUB1280: A2(J9,1 )=SS*DY 

X=XL :GOSUB1296 : GOSUB1280 :A2(J9,1 )=A2(J951 )4+SS*DY 
FORJ8=2TON: X=X+DX: GOSUB1 296 : GOSUB1280 

A2(J9>51 J=A20 59,1 )42*SS*DY :NEXTJ8 

A2(J9,1)=A2(J9,1 )*DX/3 

IFJ9=1THENNEXTJ9 

FORJ8=1TOJ9-1 

A2(J9 »J8+1 J=A2( 59,58 +l (AZ( 09 »J8 )-A2( J9-1 58) 14% J8-1 ) s NEXTIS 
T1=A2(J9,J9 J-A2( 59 ,J9-1 Ds TFSGN(T1 )*T1-TL>OTHENNEXTJ9ELSE1 266 
PRINT'Tolerance of"3TL3"“not met after five extrapolations" 
IN=A2(J9,J9):RETURN . 

FORJ7=1T06: FORJ6=1TOJ7: PRINTUSING'## .#8#" 3A2(0 57,6 J 35 

NEXTJ6: PRINT" :NEXTJ7: INPUTZ9: RETURN 

‘Simpson's Rule Sum 
Y=YU:GOSUB1 285:SS=F : Y=¥YL:GOSUB1285:SS=SS+F 

FORJ5=2TO(N/2 ): Y=Y+DY : GOSUB1 285: SS=SS+4*F : Y=Y+DY : GOSUB1285 
SS=SS+2*F ;NEXTJ5: Y=Y+DY ;: GOSUB1285: SS=SS+4*F : RETURN 

F=1 

RETURN 

XUPPER=1 

XLOWNER=0 

RETURN 

YUPPER=X 

YLOWER=0 

DY=( YU-YL JZ (N41 ):RETURN 


Figure 8.12 Integration Subroutine. 
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APPENDIX A 
DETECTION SIMULATION PROGRAM LISTING 


A complete listing of the Detection Simulation Program is as follows. 


100 CLS:PRINT"":PRINT" DETECTION SIMULATION": FORI=1T0400:NEXTI 
110 ‘Input/Initialization 

115 OPEN"DSIN"FORINPUTAS1 

120 INPUT#1 NS ,NP»,S1,S2,RH,F1:V1l=S1¥S1:V2=S2*S2 

121 DIMX(NS,5+NP ),A2(65,6),T1(3) 

125 FORIL=1TOS+NP: INPUT#1,X(1,I1):NEXTIIL 

126 IFNS=1THEN14¢0 

127 IFFI=1THEN132 

128 FORI1=2TONS: FORI2=1TO5S+4NP;: INPUT#1,X(1I1,I2):NEXTI2:NEXTI1:GOTO140 
132 FORI1=2TONS: FORI2=1TO5: INPUT#1 »X(I1,12):NEXTIZ2: IFNP=OTHEN135 

134 FORI2Z=6TOS4NP:X(I1,I12)=X(1,I12):NEXTI2 

135 NEXTI1 

140 RF=SQR(1~RH*%2) 

150 DFS$="EXP(-( (XT-XS )424(YT-YS )% 2/0 2*X( 526 )%2))" 

200 '**¥* Simulation Selection Section *** 

201 CLS:PRINT"Is the Detection function:" 

203 PRINT" 1. Oeterministic":PRINT" 2. Probabilistic" 

205 INPUT"Enter 1 or 2:"3Fl1 

210 CLS:PRINT"Are Sensor Locations:" 

212 PRINT" 1. Always At Aim Point":PRINT" 2. Distributed BYN Around Aim Point" 
214 INPUT"Enter 1 or 2:'"3F2 

215 IF( F1=2ORF2=2 JTHENF3=1:GOTO230 

220 CLS:PRINT"":PRINT"Is the Calculation:" 

222 PRINT" l=Monte Carlo Simulation":PRINT" 2=Numerical Approximation" 
224 INPUT"Enter 1 or 2"5F3 

230 TIMES="00:00:00":IF FI=1THENGOSUB300ELSEGOSUBS500 

250 GOTO200 

300 ‘Deterministic Sensor Subroutine 

305 IFFS=1THENGOSUB310ELSEGOSUB350 

306 RETURN 

310 ‘Monte Carlo of Deterministic Sensor 

315 GOSUB900 

320 PD=0: FORJI=1TONR:PRINT.241,"Repetition:"3J1:GOSUB600: FORJ2=1TONS 
323 IFF2=2THENXS=X(J2,1):YS=X(J2,2 ):GOTO325 

324 GOSUB612 

325 T1=SQR( (XS-xT )* 2+(YS-YT )%2) 

330 IFT1I<=xX(J2,6 JTHENPD=PD+1:GOTO335 

332 IFT1>=X(J2,7 JANDT1<=X(J2,8)THENPD=PD+1:GOTO335 

33¢ NEXTJ2 

335 NEXTJ1:PD=PD/NR:GOSUB950: RETURN 

350 'Numeric/Deterministic Subroutine 

355 PD=0: FORJ2=1TONS: H=X(J2,6):GOSUB1200: PD=PD+IN 

356 H=X(J2,8):GOSUB1200:PO=PO+IN 

357 H=X(J2,7):GOSUB1200:PD=PD-IN:NEXTJ2 

360 GOSUB950: RETURN 

500 'Probabilistic Detection Function 

502 CLS:PRINI"Oefault Oetection Function Is Carleton." 

503 GOSUB1300:GOSUB900 

520 PD=0: FORJI=1TONR: PRINT. 241,"Repetition:"3J1:GOSUB600 : FORJZ=1TONS 
521 IFF2=2THENXS=X(J251):YS=X(J2,2 ):GOTO523 

522 GOSUB612 

523 GOSUB1410: IFRNO(1)<=DFTHENPD=PD+1:GOTO526 

52¢ NEXTJ2 

526 NEXTJ1:PD=PD/NR 

530 GOSUB950: RETURN 

600 'x**Generate BYN RV 

602 U1=RNO(1):U2=RND( 1): TE=SQR( -2*LOG(U1 )) 
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604 
606 
612 
614 
616 
900 
905 
910 
950 
951 
952 
955 
954 
955 
956 
957 
958 
959 
960 
961 
962 
963 
965 
966 
967 
968 
970 
97¢ 
974 
976 
O77. 
1200 
1201 
1202 
1220 
1230 
1240 
1242 
1245 
1250 
1251 
1252 
Z55 
1260 
1262 
1263 
1264 
1266 
1275 
1276 
1280 
1281 
1282 
1285 
1286 
1287 
1290 
293 
1296 
i297 
1300 
1302 
1304 
1306 
1307 
1400 
1410 
1450 
1451 


XT=TEXCOS(6.2831853*U2 }: YT=RH*XT+RFXTEXSIN( 6. 2831853*U2 ) 

XT=XT¥S1: YT=YT*¥S2: RETURN 

U1=RND(1):U2=RND(1):TE=SQR(-2*LOG( U1) ) 

XS=TEXCOS( 6. 2831853*U2 ): YS=X( J2,5 )*XT4+(1-X( 52,5 )*2 * . 5XTEXSIN( 6. 2831853*U2 ) 
XS=X( S21 )4#XS*K( S253 3: YS=X( S252 J4YS*#X( 52 54): RETURN 

CLS: INPUT"Enter number of repetitions for Monte Carlo Simulation: "3NR 
RETURN 

INPUT"Hit ENTER to Continue"3Z1:RETURN 

‘Print output 

SOUND1567 ,10:SOUND1244 ,10: SOUND1046 , 10: SOUND 783 » 20 

SOUND 1046 , 10: SOUND 783 »540 

CLS:PRINT"":PRINT"Calculation Time (HH/MM/SS) = "3;TIME$: IFF3=2THEN960 
PRINT"Select Alpha for Confidence interval:" 

INPUT" Choices = .1, .05, .O1:"3AL 

IFAL=.1THENAL=1.645:GOTO960 

IFAL=.O5THENAL=1. 96:GOTO960 

IFAL=.O1THENAL=2.575:GOTO960 

GOTO954 

PRINT''x** Estimate of P(Detection) = "3:PRINTUSING” #.#####"3PD 
IFF3=1THEN965 . 

PRINT"No Confidence Interval For Numerical Approximations" 

GOTO970 

PRINT'Confidence Interval: ("3 

TE=AL*SQR( PD*¥(1-PD )/NR ): LL=PD-TE: UL=PD4TE : IFUL>1THENUL=1 

IFLL<OTHENLL=0 

PRINTUSING'##. ###8#"5LL3UL3: PRINT" )"': GOSUB9I10 

‘Confetti Approximation 

PRINT’: INPUT"Confetti approximation? O=No, 1=Yes:'"';Z9:IFZ9=OTHENRETURN 
CLS: INPUT"Enter TOTAL lethal area for ALL sensors in the pattern:";NA 
TE=NA/( 6. 283185*S1*S2 ): TE=1-( 1#SQR( 2*TE ) }XEXP(-SQR( 2*TE )) 
PRINT"**Confetti Approximation = ";TE:GOSUB910:RETURN 


‘Numberical Integration Subroutine 

D1=6 ..2831853*S1*S2*RF 

TL=.001 

CLS:PRINT'": PRINT" '!*Calculating An Integral!!":PRINT'" 


N=2:GOSUB1293 :DY=( YU-YL)/2 

FORJ9=1T06: DY=DY/2 :N=N¥2 
Y=YU:GOSUB1296 : GOSUB1280 :A2(J9,1 J=TS*DX 

Y=YL:GOSUB12 96 : GOSUB1280: A2(J9,1)=A2(J9,1)4+TS*DX 
FORJ8=2TON: Y=Y+DY : GOSUB1296 : GOSUB1280 

A2(J9,1)=A2(J9>1 )#2*TS*DKX: NEXTJS 

A2(J9,1 )=A2(J9,1 )*DY/2 

IFJ9=1LTHENNEXTJ9 

FORJ8=1TOJ9-1 

A2(J9,J841 }=A2( 59 ,J8 J4( (A209, J8 }-A2( J9-1,J8 ) 714% S8-1) )sNEXTIU8S 
T1=A2( 59539 J-A2(59,J59-1): IFSGN(T1 )*¥T1-TL>OTHENNEXTJ9ELSE1266 
PRINT"Tolerance of"3TL3"not met after five extrapolations" 
IN=A2(J9,J9): RETURN 

FORJ7=1T06: FORJ6=1TOJS7:PRINTUSING' 4. ##E" 5 A2(057,J6 )3 

NEXTJ6: PRINT" :NEXTJ7: INPUTZ9: RETURN 

REM Trapezoidal Rule Sum 
X=XU:GOSUB1286 : TS=F :X=XL: GOSUB1286: TS=TS+F 

FORJ5=2TON-1: X=X+DX : GOSUB1286 : TS=TS+F :NEXTJS5;: RETURN 

"F(X,Y) to be integrated: 

F=X* 2/V1-2*RH¥X*Y/S1/S2+Y% 2/V2 

F=( EXP(-F/2/RF%* 2 ))/D1:RETURN 

‘Limits of Integration: 

YU=X(J252 )+H: YL=X(J2,2 J-H: RETURN 

T3=SQR(H* 2-(Y-X( 52,2 ))*% 2): XU=K(J2,1)4+7T3: XL=X(J2,1)-T3 :DX=( XU-XL I/N 
RETURN 

PRINT" -Detection Fn (DF) in terms of XT, YT," 

PRINT" and Parameters XS» YS» and X(J2,6).... X(J2,5+NP):" 
PRINT’ *x* DF = "3DF$ 

PRINT"Hit ENTER For No Change or Enter New...": INPUT" DF = "3DF$ 
RETURN 

‘Tokenize DF 

BS="DF="+DF$+CHRS$(0) 

‘Tokenize/execute B$ 

BO=VARPTR(B$ }:B1=PEEK(BO+1 }+256*PEEK(BO+2 }:CALL1606,0,B1l 
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1455 CALL2499,0,63105:RETURN 
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APPENDIX B 
KALMAN FILTER PROGRAM LISTING 


A complete listing of the Kalman Filter Program 1s as follows. 


100 CLS: PRINT"**x**KALMAN FILTER***xx": PRINT" Input Data Being Read" 
110 OPEN'KALIN'"FORINPUTAS1 : ONERRORGOTO9900 

120 INPUT#1,NX,NZ: I FNX<NZTHENMD=NZELSEMD=NX 

125 DIMPH(NX>,NX) »MNONX),Q0NX NX) La ee Aaa R1(NZ»NZ) »>MUC NX) >SGONX NX ) 
126 DIMC1(MD,MD),C2(MD,MD),K(NX,NZ) 

127 DIMB1(NZ+1,NZ*2 ) 

130 FORI1=1TONX: FORIZ=1TONX: INPUT#1,PH(I1,I2):NEXTI2:NEXTI1 

132 FORI1=1TONX: INPUT#1 »MHW(I1):NEXTI1 

134 FORI1=1TONX: FORIZ=1TONX: INPUT#1,Q(I11,I2):NEXTI2:NEXTI1 

136 FORI1=1TONZ: FORIZ2=1TONX: INPUT#1,H(I1,I2):NEXTI2:NEXTI1 

138 FORI1=1TONZ: INPUT#1 ,MV(I1):NEXTI1 

140 FORII=1TONZ: FORI2=1TONZ: INPUT#1,R1(I1,I12):NEXTI2:NEXTI1 

142 FORI1=1TONX: INPUT#1 ,MU(I1):NEXTI1 

144 FORI1L=1TONX: FORI2=1TONX: INPUT#1,SG(I1,I12):NEXTI2:NEXTII1 

145 CC=0:CLS:PRINT"Initial SG As Input Check:":GOSUB532 

150 CLS:PRINT" *XXXXMEASUREMENT BLOCK 33%" 

160 PRINT"Current 4H ":GOSUB540 

162 INPUT"Enter New H ? 1=Yes, O=No:"329:IFZ9=OTHEN170O 

165 ‘Enter A New H 

167 FORI1=1TONZ: FORIZ2=1TONX 

168 PRINT"Enter Row''3I13;", Column"3I23"Of H :'"3 

169 INPUTH(I1,I2):NEXTI2:PRINT'":NEXTII1 

170 'CALC KALMAN GAIN 

171 'MULT SG H t, INTO Cl 

172 FORI1=1TONX: FORI2Z=1TONZ:C1(1I1,I2 )=0: FORI3=1TONX 

174 C1(1I1,I12 )=(SG(I1,13 )*H(I2,I3) )+C1l(I1,I2):NEXTIS:NEXTI2:NEXTI1 
180 'MULT H SG Ht », INTO C2 

182 FORI1=1TONZ: FORI2=1TONZ:C2(1I1,I2 )=0: FORI3=1TONX 

184 C2(1I1,I2)=(H(I1,13 )*C1(1I3,1I2) )+C2(011,12):NEXTI3:NEXTI2:NEXTI1 
200 'ADD R_ INTO C2 

202 FORI1=1TONZ: FORI2=1TONZ:C2(1I1,12)=C2(I1,I12)+R1(1I1,12) 

203 NEXTI2:NEXTI1 

210 ‘INVERT C2 

215 GOSUB 9800 

220 'MULT Cl C2 INTO K 

222 FORI1=1TONX: FORIZ=1TONZ:K(1I1,I2 )=0: FORI3Z=1TONZ 

224 K(I1,I2)=(C1(I1,13 )xC2(13,12) )+K(I1,12):NEXTIS:NEXTI2:NEXTI1 
250 'x**xx*XxXUPDATE MU- TO MU+xxxx 

251 "MULT H MU- INTO Cl 

252 FORI1=1TONZ:C1(I1,1 )=0: FORIZ=1TONX 

254 C1(I1,1)=(H(I1,13 )XMU(I3 ) )4+C1(I1,1):NEXTI3S:NEXTI1 

260 'ADD MV + H MU- 

262 FORI1=1TONZ:C1(1I1,1 )=C1(I1,1)+MV(I1):NEXTI1 

270 ‘INPUT A NEW MEASUREMENT 

272 CC=CC+1:CLS:PRINT"'Measurement #'"3CC3";" 

273 FORI1=1TONZ:PRINT"Enter Element'3I13;"Of Measurement: "3 

274 INPUTZ(1I1):NEXTI1 

280 ‘SUBTRACT Cl FROM Z, INTO Cl 

282 FORI1=1TONZ:C1(I1,1)=Z(1I1)-C1l(I1,1):NEXTI1 

290 'MULT K Cl INTO C2 

292 FORI1=1TONX:C2(1I1,1)=0: FORIZ=1TONZ 

294 C2(I1,1)=(K(I1,13 )*C1(1I3,1))4+C2(I1,1):NEXTIS:NEXTI1 

300 ‘ADD C2 + MU- TO UPDATE TO MU+ 

302 FORI1=1TONX:MU(I1 )=C2(1I1,1)+MU(I1):NEXTI1 

320 'MULT K H & SUBTR FROM I , PUT IN Cl 

322 FORI1=1TONX: FORI2=1TONX:C1( 11,12 )=0: FORI3=1TONZ 

324 C1(I1,I2)=(K(I1,13 )¥H(I3,12) )4+C1l(I1,I12):NEXTI3:C1l(I1,I2 )=-Cl( 11,12) 
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NEXTI2:NEXTI1 

FORI1=1TONX:C1(I1,11)=14+C1(I1,I1):NEXTI1 

"MULT LAST RESULT BY SG , INTO C2 

FORI1=1TONX: FORI2=1TONX:C2(I1,I2)=0: FORIZ=1TONX 
C2(I1,I23=(C1(I1,13 )¥SG(1I3,I12) 34C2(I1,I2):NEXTI3:NEXTI2:NEXTI1 
"PUT C2 INTO SG 

FORI1L=1TONX: FORI2=1TONX:SG(I1,I2)=C2(I1,I2):NEXTI2:NEXTI1 
CLS:PRINT'Kalman Gain, Kfi,j) After" 

PRINT"Measurement #''5CC:GOSUB510 

CLS: PRINT"“Estimate Of System State, MU(i)+ After" 
PRINT"Measurement #''3;CC:GOSUB520 

CLS:PRINT"Estimate Of Covar, SG(1,j)}+ After" 

PRINT'Measurement #''3CC :GOSUB530 

CLS: PRINT" **x*****XMOVEMENT BLOCK xxx" 

"Update MU(CC)+ to MU(CC4+1)- 

"MULT PH MU , PUT IN Cl 

FORI1L=1TONX:C1(I1,1 )=0: FORI3=1TONX 

C1(I1,1 )=(PH(I1,13 }x*MU(I3) )4C1(I1,1):NEXTI3:NEXTI1 

"ADD Cl+ MW », INTO WMV 

FORIL=1TONX:MU(TI1 3=C1(I1,1)4MW(I1):NEXTI1 

"xXUPDATE SG ** 

‘MULT T SG , INTO Cl 

FORI1=1TONX: FORIZ=1TONXK:C1(I1,I2 )=0: FORI3=1TONX 
C1(I1,I12)=(PH(I1,I3 )*SG(I3,I2))4C1(I1,I12):NEXTI3:NEXTI2:NEXTI1 
"MULT Cl PH t, INTO C2 

FORI1=1TONX: FORI 2=1TONX:C2( 11,12 }=0: FORIZ=1TONX 
C2(I1,I2)=(C1(I1,13 )*PH(I2,I13) }4C2(I1,I12):NEXTI3S:NEXTI2:NEXTI1 
"ADD C2 +Q= SG 


FORI1=1TONA: FORIZ=1TONK: SG(I1,12 }=C2(I1,12 }+Q(I1,12}:NEXTI2Z:NEXTI1 


PRINT'Estimate Of System State, MU(I)-" 
PRINT''Before Measurement #"3CC+1:GOSUB520 
CLS:PRINT“Estimate Of Covar, SG(I,J)}- Before" 
PRINT''Measurement #''3;CC+1:GOSUB530 
GOTO160 
‘PRINTING SUBROUTINES 
"PRINT KALMAN GAIN, K 
FORI1=1TONX: FORIZ=1TONZ: PRINTUSING" #88884. #4" 5K0I1,I12)3:NEXTI2 
PRINT" sNEXTIL: INPUT"Hit ENTER To Continue: "329: RETURN 
‘PRINT MU 
FORIL=1TONX: PRINTUSING'S#8##8#. 4H" 5MUCI1)5:NEXTI1: PRINT" 
INPUT"Hit ENTER To Continue: ''3Z9: RETURN 
"PRINT COVAR MATRIX, SG 
FORIL=1TONX: FORI2Z=1TONX: PRINTUSING#H8## . #4" 3S6(I1,I12)3:NEXTI2 
PRINT" :NEXTIL: INPUT"Hit ENTER To Continue: "3Z9:RETURN 
‘PRINT 4H 
FORI1=1TONZ: FORI 2=1TONX: PRINTUSING'###8#8# .#8'"53H(I1,I2)3:NEXTI2 
PRINT’ :NEXTI1: RETURN 
PRINT" C2 MATRIX:" 
FORI1=1TOA: FORI2=1TOB: PRINTUSING###88#. F#"5C20I11,I2)5 :NEXTI2Z 
PRINT" :NEXTIL: INPUT"Hit ENTER To Continue: ''3Z9: RETURN 
"INVERT C2 
FORIL=1TONZ: FORI2=1TONZ:B1(I1,I2  }=C2(I1,I2):NEXTI2:NEXTI1 
FORI1=NZ+1TO2¥NZ: FORI2=1TONZ 
IFIL=I24+NZTHENB1(I2,I1)J=lLELSEB1(I2,I1}=0 
NEXTI2:NEXTI1 
FORI1=1TONZ 
ML=1/B1(I1,I1):FORI3Z=1TO2*¥NZ:B1(I1,1I3 }=B1(I1,13 J*ML:NEXTI3 
IFIL=NZTHEN9865 
FORI2=I11+1TONZ: IFB1(1I2;,I1 )J)=OTHEN9860 
ML=-B1(I2,I1)} 
FORI3Z=11TO2*NZ:B1(1I2,13 J=B1(12,1I3 }+(ML¥B1(I1,I3)):NEXTI3 
NEXTI2:NEXTI1 
FORI1=NZTO2STEP~1 
FORI2=I1-1TO1LSTEP-1:IFB1(I2,I1 J=OTHEN9885 
ML=-B1(1I2,I1) 
FORI3=1TO2*NZ:B1(I2,1I3 J=B1(1I2,I3 J4(ML¥B1(I1,I3)):NEXTI3 
NEXT12:NEXTI1 
FORI1=1TONZ: FORI2=1TONZ 
C2(I2,I1}J=B1l(IZ2,I14NZ):NEXTI2Z2:NEXTI1 
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9897 MI=1:RETURN 
9900 IFERL>9700ANDERR=1LITHENPRINT'!!!'ERROR: C2 Is Not Invertable!!!":END 
9905 PRINT"Error Code"3ERR3;"In Line" 3;ERL:END 
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APPENDIX C 
LANCHESTER SIMULATION PROGRAM LISTING 


A complete listing of the Lanchester Simulation Program is as follows. 


100 ‘LANCHESTER TIME STEP MODEL 

120 MAXFILES=2:OPEN"LANIN” FORINPUTAS1 

121 OPEN'LANOUT"'FOROUTPUTAS2: INPUT#1 ,NP,NA,ND 

122 IFNA>NDOTHENMD=NAELSEMD=ND 

124 IFMD<STHENCLS:SF=1 . 

130 DIMAA(NA,ND),BB(ND»NA),AT(NA),DT(OND),AR(NA),DROND) 

131 DIMQA(2,NA),QD(ND),AB(2,NA),0B(2,ND) ,SAC(NA) »SD(ND) ,OA( NA) ,OD( ND) 
132 ‘Enter Initial Quantities of Wpns, Break Points And Wpn Types. 
134 FORI2Z=1TONA: INPUT#1,QA(2,12):SA(I2 J=QA(2,1I12): OAC I2 J=127:NEXTI2 
135 FORI2=1TOND: INPUTS1,QD(1I2):SD(I2 J=QD(I2 ):O00( I2)=238:NEXTI2 

136 FORI2=1TONA: INPUT#1 ,AB(1,I12):AB( 2,12 J=AB( 1,12 )¥QA( 2,12): NEXTI2 
137 FORI2=1TOND: INPUT#S1,D8(1,12):08(2,I12 )=DB( 1,12 )*QD(1I2):NEXTI2 

138 FORI2=1TONA: INPUT#1,AT(I2 ):NEXTI2: FORIZ=1TGND: INPUT#1,DT(I2 ):NEXTI2 
140 TM=0: IFSF=1THENGOSUB600 

143 FORIL=1TONP:PRINT#2,"STARTING PHASE"3;I1 

145 ‘Enter Time Spent In Phase Il and # of Intervals 

146 INPUT#1,TT»NI:OT=TT/NI 

150 ‘Enter Replacement Rates And Attrition Coefficient Matrices 

152 FORIZ=1TONA: INPUT#1,AR(I2):NEXTI2: FORIZ=1TOND: INPUT#1,DOR(1I2 ):NEXTI2 
154 FORI2=1TONA: FORI3=1TOND: INPUT#1,AA(I2,I3) 

156 NEXTI3:NEXTI2 

157 FORI2=1TOND: FORI3=1TONA: INPUT#1,BB(1I2,I3) 

159 NEXTI3:NEXTI2 

200 ‘Fight Phase Il. 

202 FORIZ2=1TONI: TM=TM+DT:PRINT.291,"Phase:"3113", Increment"31I23;"out 
of";NI 

210 ‘Fight Time Increment DT. 

220 ‘Update number of attackers 

222 FORIZ=1TONA: QA(1,I3 J=QA( 2,13 ):NEXTI3: FORIS=1TONA: FORIG=1TOND 

223 ‘IFDOT(IG¢ J=1LTHENQA( 2,13 J=QA( 2,13 )-AA( 13,14 )¥QD( 14 )¥QA( 2,13 )¥DT: GOTO226 
224 ‘QA 2,13 )=QA( 2,13 )-AA(I3,14¢ )*¥QD( 14 J¥DT 

225 QA( 2,13 J=QA( 2,13 )-AA(I3 514 )*(QA( 2,13 )/QD( 14) )“DT( 14 )*QD( 14 )*DT 
226 NEXTI4:QA( 2,13 J=QA( 2,13 J+AR(I3 )¥DT: IFSF=1THENGOSUB650 

227 NEXTIS 

230 ‘Update number of defenders 

232 FORI3S=1TOND: FORIG=1TGNA 

233 ‘IFAT( IG )J=1THENQD(I3 )=QD( 13 )-BB( 13,14 )*¥QD( 13 )¥QA(1,14¢)*DT: GOTO236 
234 *QD( 13 )=QD0(I3 )-BB( 13,14 )*QA( 1,14 )*DT 

235 QD( 13 )=QD( 13 )-BB( 13,14 )*( QD( 13 )/QA( 1,14) )*AT(14 )¥QA(1,14)*DT 
236 NEXTI4: QD( 13 J=QD( I3 )+DR( 13 )¥OT: IFSF=1THENGOSUB660 

237 MEXTIS 

240 GOSUB300:NEXTI2 

242 IFI1L=NPTHENGOSUB350:CLS:PRINT"Output is in file LANOUT.DO.":END 
245 PRINT#2,"Status After Phase'";I1:GOSUB361:NEXTI1 

300 ‘Check Whether Breakpoint is reached. 

320 TF=0: FORIS=1TONA: IFQA( 2,13 )>AB( 2,13 JTHEN3S25 

322 TF=1:PRINT#2,"“Attacker Wpn";133;"Is Below Breakpoint" 

323 PRINT#2," Bp ="3:PRINT#2 USING #8845." 5AB(2,13 )3 

324 PRINT#2," Current Level ="3:PRINT#2 ,USING"##8S .##"53QA( 2,13) 

325 NEXTI3 

335 FORI3=1TOND: IFQD( 13 )>DB( 2,13 )THEN3S40 

336 TF=1:PRINT#2,"Defender Wpon";I33;"Is Below BreakKpoint" 

337 PRINT#2," Bp ="5:PRINT#2,USING'#hHe. #2"50B(2,13 )3 

338 PRINT#2>," Current Level ="3;:PRINT#2,USING'SHH84 . $3" 5;QD(13 ) 

340 NEXTI3: IFTF=OTHENRETURN 

350 PRINT#2,"":PRINT#2Z,"": PRINT#2>"SUMMARY AT END OF BATTLE" 

351 PRINT#2,"":PRINT#2,"Time Elapsed During Battle ="$ 


; nO? 


352 
355 
361 
363 
364 
366 
367 
368 
600 
610 
620 
623 
625 
627 
628 
630 
632 
633 
635 
650 
653 
655 
656 
657 
658 
659 
660 
663 
665 
666 
667 
668 
669 


PRINT#2 ,USING"#### . #4" 3TM: PRINT#Z,""': GOSUB361 
CLS:PRINT" Output is in file LANOUT.DO";END 

PRINT#2," Att Wen Breakpoint Current Level" 
FORI3=1TONA: PRINT#2Z USING ######8" 5135 

PRINTH2Z USING "H#H#HHHHEH . HH" 5AB( 2,13 )3;QA( 2,13) :NEXTI3:PRINT#Z,"" 
PRINT#2,"" Def Wpn Breakpoint Current Level" 
FORI3=1TOND: PRINT#2 ,USING' #H#HHHE" 51335 

PRINT#2 ,USING" ####HHHHEE . HH" 5DB( 2,13 )5;QD(13 ):NEXTI3: PRINT#Z,“"*: RETURN 
‘Set up output screen 

PRINT"Wpn # Attacker Defender" 
FORI1=1TOMD:PRINTUSING'##" 311 

TP=2+1I11%8 

IFI1L>NATHEN630 

LINE(18,TP )-(119,TP+4),1,B:BP=18+INT(100%AB(1,I1)) 
LINE(BP-1,TP+1 )-(BP,TP+3),1,B 

IFI1L>NDTHEN635 

LINE(138,TP )-(239,TP+4),1,B:BP=138+INT(100*DB(1,I1)) 
LINE(BP-1,TP+1 )-(BP,TP+3),1,B 

NEXTI1: RETURN 

‘Update screen output of attackers 

TP=3+1I3%*8 

LINE(OA(I3),TP)-(OA(I3),TP+2),0 

OA(I3 J=18+INT( 100*QA( 2,13 )/SA(I3)) 

IFOA(I3 )>LISTHENOA( 13 J=118: PRINT. (13*40+2),°"'*'':GOTO659 
PRINT .I13*40+2," " 

LINE( OA(I3),TP )-(OA(I3 ),TP+3 ),1: RETURN 

‘Update screen output of defenders 

TP=3+13*8 

LINE(OD(I3),TP )J-(OD(I3),TP+2),0 

OD( 13 )=138+INT(100*QD(1I3)/SD(I3)) 

IFOD( 13 )>238THENOA( I3 )=238: PRINT .I13*40+22,"%":GOTO669 
PRINT .13*40+22," " 

LINE(OD(I3),TP)-(OD(1I3),TP+3 ),1:RETURN 
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APPENDIX D 
GEOMETRIC PROGRAMMING PROGRAM LISTING 


A complete listing of the Geometric Programming Program is as follows. 


100 ‘Geometric Programming Program 

110 OPEN"GEOIN"FORINPUTAS1 

120 INPUT#1>NT »NV:K9=NT 

122 IFNT-NV<>1THENPRINT"***ERROR: Degree of Difficulty <> O”":END 
130 INPUT#1,NC:DIMNT(NC ) : 

140 MN=0:FORI1L=OTONC: INPUT#1 NTC I1): IFNTCI1)>MNTHENMN=NT(I1):NEXTI1 
143 DIMCT(3,NC,MN),LM(NC ),B1(NT #1 ,NT¥2),B2(NT ) »B3CNT NV) 

145 FORI1=OTONC: FORIZ=1TONT(1I1):INPUT#1,CT(1,I1,12):NEXTI2:NEXTI1 
150 FORI1=1TONT(0):B1(1,I1)=1:NEXTI1 

155 FORI1=NT(0)+1TONT:B1(1,I1 )=O:NEXTI1 

160 FORIL=2TONT: FORIZ=1TONT : INPUT#1,B1(I1,I2):B3(1I2,I1-1)=Bl(I1,1I2) 
162 NEXTI2:NEXTI1 

170 PRINT": PRINT" **XCOMPUTING DELTA'S*x" 

172 B2(1)=1:FORI1=2TONT :B2(1I1 )=0:NEXTI1 

180 GOSUB9800 

200 CLS:I1l=1:FORI2=OTONC: FORI3=1TONT(1I2):CT(2,1I12,13 J=B1l(1I1,1) 

203 PRINT"DELTA('3I23","3I335") = "5 

204 PRINTUSING' #Hh88. ##HH"5CT(2,12,13 ): IL=I1+1:IFIL>STHENGOSUB6O00 
205 NEATI3:NEXTI2: GOSUB600:CLS 

210 PRINT": PRINT" **COMPUTING OPT FN VALUEXx*" 

212 FORI1=OTONC: LM(I1)=0:FORIZ2=1TONT(11):LM(I1 )=LM(11)+CT(2,11,12) 
214 NEXTI2:NEXTI1 

220 FS=1:FORI1=OTONC 

222 FORIZ=1TONT(1I1):FS=FS*¥(CT(1,11,12)/CT(2,11,12) )*CT(2,11,12) 
224 NEXTI2: FS=FS*¥(LM(I1)“LM(1I1)):NEXTI1 

229 PRINT" :PRINT"FX ="3:PRINTUSING' #e8#.###"5FS:GOSUB600:CLS 

230 ‘Compute optimal x(n) 

232 K9=K9-1 

234% FORIL=1TOK9: FORIZ=1TOK9:B1(1I1,I2)=B3(I1,I2):NEXTIZ2:NEXTI1 

236 CC=1: FORIL=OTONC: FORIZ=1TONT(I1) 

237 CT(3,I1,12)=(CT(2,11,I12)/CT(1,11,I12)/LM(I1)) 

238 IFI1=OTHENCT(3,I1,12)=CT(3,I1,12)J*FS 

239 B2(CC )=LOG(CT(3,I1,12)):CC=CC+1:NEXTIZ2:NEXTI1 

242 PRINT"P(m,t)* = opt. value of term t, constr. m, divided by its coefficient." 
2464 FORI1L=OTONC: FORIZ=1TONT(I1):PRINT"P("5I15","5125" )* ="'5 

2466 PRINTUSING'#H#83 . #888"5CT(3,I11,12):NEXTI2: GOSUB600:NEXTI1:CLS 
250 PRINT" :PRINT"** Computing Opt Values Of X{(n) **" 

260 GOSUB9800 

270 CLS: FORI1=1TOK9:PRINT"X¥("3I15") = "3 

272 PRINTUSING'RSRSHF .FHEHH" 5 EXP(B1(I1,1)) 

273 IFI1L>STHENGOSUB6O00 

275 NEXTI1:GOSUB600:END 

600 INPUT"x* Hit ENTER To Continue: "3Z9:RETURN 

9800 'Simultaneous Linear Equation Subroutine: Ax=b 

9815 ‘Invert Matrix A 

9820 FORK7=K9+1TO2*K9: FORK8=1TOK9 

9822 IFK7=K8+K 9THENB1(K8>K7 J=LELSEB1(K8 »K7 )=0 

9825 NEXTK8:NEXTK7 

9830 FORK7=1TOK9 

9835 IFB1(K7,K7 )*SGN(B1(K7,K7) )<1E-8THENGOSUB9I910 

9840 K2=1/B1(K7,K7):FORK6=1TO2*K9:B1(K7,K6 )=Bl1(K7,K6 J¥K2:NEXTK6 
9842 IFK7=K9THEN9865 

98465 FORK8=K7+1TOK9: IFB1(K8,K7 J)=OTHEN9860 

9850 K2=-B1(K8,K7) 

9855 FORK6=K7T02*K9:B1(K8 >K6 )=B1(K8 >K6 J +(K2*B1(K7,K6 ) ) SNEXTK6 
9860 NEXTK8:NEXTK7 

93865 FORK7=K9TO2ZSTEP-1 
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9870 
9875 
9880 
9885 
9890 
9894 
9896 
9900 
9903 
9905 
9910 
9915 
9920 
9930 
9940 


FORK8=K7-1TOLSTEP-1:IFB1(K8,K7 )=OTHEN9885 

K2=-B1(K8>5K7 ) 

FORK6=1T02*K9:B1(K8,K6 J)=B1(K8,K6 )+(K2*B1(K7 K6 ) )s NEXTK6 

NEXTK8:NEXTK7 

‘Mult A Inverse by b 
FORK7=1TOK9:B1(K7,1 )=0: FORK8=1TOK9:B1(K7,>1 )=B1(K7,1 )+B1(K7,K8+K9 )¥B2(K8 ) 
NEXTK8:NEXTK7: RETURN 

‘Error Routine 

IFERL>9700ANDERR=1L1THENPRINT"! "ERROR: Matrix Is Not Invertable!!!";:END 
PRINT"Error Code";ERR3"In Line";ERL:END 

‘SWITCH RONS 

FORK5=K7+1TOK9: IFB1(K5>,K7 )*SGN(B1(K5 ,K7 ) )<1E-8THEN9940 

FORKG=1TOK9*2: K3=B1(K7>,K4):B1(K7 >K% J=B1(K5 >KG ) 
B1(K5,K4 )=K3 :NEXTK4: RETURN 

NEXTKS:PRINT"“Error: Matrix Not Invertable":END 
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100 
105 
107 
300 
305 
510 
501 
504 
505 
506 
508 
509 
510 
512 
513 
514 
Bits 
516 
Say 
518 
520 
524 
526 
530 
531i 
532 
540 
548 
549 
550 
560 
570 
700 
702 
800 
810 
Biz 
815 
820 
825 
1000 
1005 
1007 
1008 
1010 
1020 
1021 
1022 
1023 
1025 
1030 
1040 
1060 
1080 
1081 
1082 
1090 
1900 


APPENDIX E 
MATRIX ALGEBRA PROGRAM LISTING. 


A complete listing of the Matrix Algebra Program is as follows. 


CLS:PRINT““:PRINT" **¥* MATRIX ALGEBRA PROGRAM *xx"';PRINT'' 
PRINT'IS INPUT MATRIX, 'MATIN.DO’ IN RAM?": INPUT" O=NO, 1=YES"3FF 
IFFF=1THENOPEN'MATIN"FORINPUTAS1 

PRINT" *x*Enter The Single Largest Dimension of" 

INPUT"The Largest Matrix To Be Processed: "3K 
DIMA1(3>,K,K),B1(K4+1,K¥2 ) »R(4),C( 4) ,DET( 2 ):MI=1:0F=1:SF=0 

CLS: EF=O0;:PRINT"*x***XMATRIX ALGEBRA PROGRAM MENUxxxx"" 


PRINT" 1. Enter Starting Left Side Matrix" 

PRINT" 2. Matrix Inversion" 

PRINT" 3. Matrix Addition":PRINT" G. Matrix Multiplication" 
PRINT" 5. Simultaneous Linear Equations" 

PRINT" 6. Print Current Answer Matrix'":PRINT" 7. Other Options" 
INPUT" **Enter Number:"3CH 


IFCH=1THENMI=1:29=0 : GOSUB7006 

I FCH=2THENMI=1 : GOSUB2000 

IFCH=3THENGOSUB3000 

IFCH=4THENGOSUB5000 

IFCH=5THENGOSUBSOO0O0 

IFCH=6THENGOSUB6000 

IFCH<>7THENGOTOSO1 

CLS:PRINT"***MORE CHOICES:":PRINT" 1. Determinant" 

PRINT" 2. Matrix Integer Exponentiation" 

PRINT" 3. Store Current Matrix" 

PRINT" 4. Retrieve Stored Matrix" 

PRINT" 5. Scalar Multiplication" 

PRINT" 6. Other Options": INPUT"x**Enter Number: "3 CH 

IFCH=1THENMI=1:GOSUB1000 

IFCH=2T HENMI=1:GOSUB7600 

IFCH=3THENGOSUB8000 

IFCH=4THENMI= = GOSUB8200 

IFCH=5THENMI =1 : GOSUB5100 

GOTOS501 

*PAUSE CONTROL 

INPUT"** Hit Enter To Continue"32Z9:RETURN 
‘INTERMEDIATE MODIFICATIONS 

PRINT'**Modify The 2nd Matrix?" 

INPUT" O=No, 1=Invert, 2=Scalar Multiply: "329 

IFZ9=OTHENRETURN 

MI=2:IFZ9=1THENGOSUB2000ELSEGOSUB5100 

GOTO810 
‘CALC DETERMINANT 

IFR(MI J=C(MI JTHENI1O1O 

PRINT"ERROR: Number of rows/columns not equal:" 
PRINT" MATRIX IS NOT INVERTABLE!":GOSUB700:EF=1:RETURN 
FORI1=1TOR(MI }: FORIZ2=1TOC( MI ):B1(I1,I2 J=AL(MI ,I1,I2 ):NEXTI2:NEXTI1 
DET(MI J=1: FORII=1TOR(MI ) 

IFB1(I1,I1 J*¥SGN(B1(I1,I1) J)<1E~LOTHENGOSUB190O0ELSE1023 
IFEF=1THEN1008 

DET(MI J=DET(MI )¥B1(I1,I1):IFI1L=R(MI JTHEN1O80 
FORI3Z=1TOC(MI ):B1(I1,I3 J=B1(I1,I3 )/B1(I1,I1):NEXTI3 
FORI2=I1+1TOR(MI }: IFB1I(I2,I1 J=OTHEN1060 
FORI3Z=I1TOC(MI J:B1(1I2,13 }=B1(I2,13 )-(B1l(I2,I1)*B1(I1,13) ):NEXTIS 
NEXTI2:NEXTI1 

IFOF<>1THENRE TURN 

PRINT"**Det. Of Matrix "3MI;" Is: "$5 
PRINTUSING'S###8# . SHEER" 5 DETCMI ): GOSUB700 

RETURN 
‘SWITCH ROWS 


Til 


1910 FORJ=I141TOR(MI ): IFB1(J,I1)*SGN(B1(J,I1) )<1E-1LOTHEN194¢0 

1920 FORJ1=1TOC(MI )¥*2: TE=B1(I1,J1):B1(I1,J1)=Bl1(J,J1) 

1930 B1(J,J1)=TE:NEXTJ1:GO0TO1950 

1940 NEXTJ:EF=1: RETURN 

1950 DET(MI J=-DET(MI ): RETURN 

2000 ‘MATRIX INVERSION 

2010 OF=0:GOSUB1000: IFDET(MI J)¥SGN( DET(MI ) )>1LE-1LOOREF=1THEN2017 

2015 PRINT"XERROR: Determinant=0. MATRIX NOT INVERTABLE!":GOSUB700:EF=1 
2017 IFEF=1LTHENRETURN 

2020 FORIL=1TOR(MI ): FORIZ=1TOC(MI ):B1(I1,I2 J=A1L(MI,I1,I12):NEXTI2:NEXTI1 
2030 FORIL=C(MI )41TO2*C( MI): FORIZ=1TOR(MI ) 

2032 IFI1L=I12+R(MI JTHENB1(I2,I1)J=LELSEB1(I2,I1)=0 

2035 NEXTI2:NEXTI1 

2040 FORI1=1TOC(MI) 

20465 IFB1(I1,I11)*SGN(B1(I1,I1) )<1LE-LOTHENGOSUB19SOOELSE2055 

2046 IFEF=1THEN2Z015 

2055 ML=1/B1(1I1,11):FORIZ=1TO2*C(MI ):B1(I1,13 )=B1(I1,I3)*ML:NEXTI3 
2057 IFI1=C(MI JTHENZ080 

2060 FORI2Z=I1+1TOR(MI ): IFB1(I2,1I1)=OTHEN2075 

2065 ML=-B1(1I2,I1) : 

2070 FORI3Z=I1TO2*C(MI ):B1(I2,I3 )=B1l(1I2,1I13 J4(ML*¥B1(I1,I3)):NEXTI3Z 
2075 NEXTI2:NEXTI1 

2080 FORI1=C(MI JTO2STEP-1 

2100 FORI2=I1-1TO1STEP-1:IFB1(1I2,I1 )=OTHEN2130 

2110 ML=-B1(I2,I1) 

2120 FORI3=1TO2*C(MI ):B1(1I2,13 J=B1l(I2,13 J+(ML*¥B1(I1,I3 )):NEXTIS 
2130 NEXTI2:NEXTI1 

2140 FORI1=1TOC(MI ): FORI2Z2=1TOR(MI ) 

2145 A1L(MI,I2,11)=B1(I2,I14C(MI) ):NEXTI2:NEXTI1 

2190 OF=1:RETURN 

3000 ‘MATRIX ADDITION 

3010 MF=2:MI=2:GOSUB7000 :GOSUB800: IFEF=1THENRETURN 

3015 FORI1=1TOR( 1): FORI2Z=1TOC(1):A1(1,1I1,I2)=A1(1,I1,12)+A1(2,I1,12) 
3020 NEXTI2:NEXTI1:GOSUB6000:MF=0: RETURN 

GO00 ‘SIMULTANEOUS LINEAR EQUATIONS 

4010 CLS: PRINT"**Solves Ax=b. Choices:":PRINT" 1. Enter b Vector." 
G012 PRINT" 2. Change An Element In Matrix A." 

4013 PRINT" 3. Solve Current Ax=b." 

4014 PRINT" 4. Return.": INPUT" * Select A Number: "3CC 

4020 IFCC=1GOTO4¢040 

GO022 IFCC=2GOTOG050 

G024 IFCC=3GOTOG060 

G026 IFCC=G4THEN RETURN 

4035 GOTO 4000 

G0G0 MI=2:R(2 J=C(1):C( 2 J=1:GOSUB7040: GOTOG000 

4050 INPUT'**xRow, Column Of Matrix A To Be Changed: '3RD,CD 

GO52 PRINT" - Enter Row"3RD3", Column";CD3":''3 sINPUTA1(1,RD,CD):GOTOG000 
4060 MI=1:SF=1:0OF=0:GOSUB8000 : GOSUB2000 

G06G IFEF=O0THENGO70 

G065 PRINT"*Solution Not Uniquely Determinable'":GOSUB700: RETURN 
4070 GOSUB5000:CC=0: FORII=1TOR( 2 ):CC=CC+#1:PRINT"x("3I13;") = "3 

G075 PRINTUSING' HHH. HHH" 581(01I1,1): IFCC>6éTHENGOSUB700:CC=0 

4080 NEXTI1:GOSUB700:SF=0:GOSUB8200 : GOTO4¢000 

5000 ‘MATRIX MULT 

5010 MF=1:IF SF=1THEN5020 

5015 MI=2:GOSUB7000 :GOSUB800: IFEF=1THENRETURN 

5020 R(G4)=R(1):C(4 J=C( 2): FORIL=1TOR(G ): FORIZ=1T0C(4 ):B1(I1,I2)=0 
5022 FORI3=1T0C(1):B1(I1,I12)=A1(1,1I11,13 )¥A1(2,1I13,I2)4+B1(I1,I2) 

5024 NEXTI3:NEXTIZ2:NEXTI1:MF=0 

5050 IF SF=OTHENGOSUB7500 :GOSUB6000 

5060 RETURN 

5100 ‘SCALER MULT 

5110 INPUT"Enter Scalar Multiplier:"3SM:FORI1=1TOR(MI ): FORI2=1TOC(MI ) 
5115 A1(MI,I1,I2 J=AlL(MI,I1,I2)*SM:NEXTI2:NEXTI1: RETURN 

6000 ‘PRINT OUTPUT MATRIX 

6010 PRINT" ** Current Answer Matrix:":CC=0;:FORI1L=1TOR(1):CC=CC+1 
6012 FORI2=1TOC(1):PRINTUSING"#HHR. HAHAH" 5A1°01,I11,12)3:NEXTIZ2: PRINT" 
6050 IFCC=STHENGOSUB700:CC=0 

6060 NEXTI1:GOSUB700: RETURN 


112 


7000 
7001 
7002 
7003 
7004 
7006 
7007 
7008 
7009 
7010 
7011 
7012 
7013 
7014 
7015 
7017 
7018 
7020 
7021 
7022 
702% 
7026 
7030 
7031 
7032 
7034 
7035 
7036 
7037 
7040 
7042 
7044 
7050 
7052 
7500 
7510 
7512 
7600 
7610 
7620 
7622 
7630 
8000 
8010 
8012 
8200 
8210 
8212 


"MATRIX INPUT 

CLS:PRINT'':PRINT"Will This Matrix Be On:" 

INPUT" O=Left, 1=Right"3Z9: IFZ9=1THEN7006 

R( 2 J=RO1):C0 2 J=C( 1): FORIL=1LTOR( 1): FORI2Z=1TOC( 1) 
A1(2,I1,I2 J=A1(1,I1,I2):NEXTIZ2:NEXTI1:MI=1 - 

CLS: IFMI=2THEN7008 

PRINT"**Choices For Left Hand Matrix':GOTO7009 
PRINT"**Choices For Risht Hand Matrix:" 

PRINT" 1. Enter Matrix From Keyboard" 

PRINT" 2. Retrieve Stored Matrix": IFFF<>1THEN7012 
IFEOF(1)=OTHENPRINT" 3. Enter Matrix From MATIN.DO" 
INPUT"*x*Enter A Number: "3Z9:IFZ9=1THEN7015 
IFZ9=3THEN7018 

R(MI J=R(3 9: COMI J=C(3):GOTO7020 

PRINT" **Enter The Rows, Columns" 

INPUT"In The Next Matrix: "3R(MI),C(MI ):GOTO7020 
INPUT#1,R(MI ),C(MI ) 

IF MF<>1THEN7030 

IFR( 2 J=C( 1 JTHEN7030 

PRINT"**ERROR: Columns in LEFT MATRIX ="3C(1) 

PRINT" Rows In Risht Matrix ="3R(2) 

PRINT"These Must Be Equal For Matrix Mult!!";:GOSUB700:EF=1:GOTO7006 
IF MF<>2THEN7035 

IF(R(1J=R( 2 VANDC(1)=C(2))THEN 7035 

PRINT" **XERROR:Dimensions For Both Input" 
PRINT"Matrices Must Be Equal!!":GOSUB700:EF=1:GOTO7006 
IFZ9=2THENGOSUB8200: RETURN 

IFZ9=3THEN7050 

PRINT" **Fill Matrix Row By Row:":PRINT"" 
FORI1=1LTOR(MI J): FORIZ=1TOC(MI ) 

PRINT" =-Enter Row"3I13"And Column''3123"':''3 
INPUTAIL(MI,I1,I2):NEXTI2:PRINT’":NEXTI1: RETURN 
FORII=1TOR(MI J): FORI2=1TOC( MI): INPUT#1,A1L(MI,I1,1I2) 
NEXTI2:NEXTI1: RETURN 

"MAKE ANSWER MATRIX THE FIRST MATRIX FOR THE NEXT OPERATION 
R(1)J=R(4):C( 1 J=C(4):3 FORILELTOR(1 J: FORI2Z=1TOC(1) 
Al(1,I1,1I12)=Bl(I1,I2):NEXTI2:NEXTI1: RETURN 

"MATRIX INTEGER EXPONENTIATION 

CLS:PRINT"': INPUT"*xEnter Positive Integer Exponent: '3XP 
R(2)J=R(1):C( 2 J=C(1):FORIL=1TOR( 1 ):FORIZ=1TOC( 2 ) 
Al(2,I1,I2)J=A1(1,I11,12 ):NEXTI2:NEXTI1 

SF=1:FOREX=2TOXP : GOSUB5020 : GOSUB 7500 :NEXTEX: GOSUB6000 :SF=0: RETURN 
"STORE Al(1,>) 
R(3J=R(1):C(3 J=CC 1): FORIL=1TOR( 3): FORIZ=1TOC( 3) 
A1l(3,I1,I2)=A1l(1,I11,I2 ):NEXTI2:NEXTI1: RETURN 

"RETRIEVE THE STORED MATRIX 

R(MI J=R03):CCMI J=C( 3): FORIL=1LTOR(MI J): FORIZ=1TOC(MI ) 
A1(MI,I1,I2J=A1(3,I1,I2 ):NEXTI2:NEXTI1: RETURN 
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APPENDIX F 
NUMERICAL INTEGRATION PROGRAM LISTING 


A complete listing of the Numerical Integration Program js as follows. 


1200 'Numberical Double Integration:Steven H. Cary:24¢ Apr 86 

1201 DIMA2(6,6):TL=.001 

1205 CLS:PRINT'":PRINT" **% Double Integration *x" 

1206 PRINT" Romberg Algorithm" 

1210 PRINT"O=Edit Function To Be Integrated." 

1211 PRINT"1=Edit Limits Of Integration." 

1213 PRINT"2=Edit Tolerance; Current Tol.="3TL 

1215 PRINT"3=Calculate Integral. ":INPUT"Enter 0, 1, 25 or 3:"329 
1216 IFZ9=OTHENEDIT1285-1288 

1217 IFZ9=1THENEDIT1291-1298 

1218 IFZ9=2THENINPUT"Tolerance = "3TL:GOTO1205 

1220 CLS:PRINT"":PRINT" 'IBe Patient!!";:PRINT"" 

1230 N=2:GOSUB1293 : DX=( XU-XL )/72 

1240 FORJ9=1T06 : DX=DX/2 : N=N¥2 

-124¢2 X=XU:GOSUB1296 : GOSUB1280:A2(J9,1)=SS*DY 

1245 X=XL:GOSUB1296: GOSUB1 280: A2(J9>51 J=A2(J9,1)4SS*DY 

1250 FORJ8=2TON: X=X+DX:GOSUB1296 :GOSUB1280 

1251 A2(J9>1)=A2(J9,1 )+2*%SS*DY :sNEXTIJ8 

1252 A2(J9s1)=A2(J9>1 )*DX/3 

1255 IFJ9=1THENNEXTJ9 

1260 FORJ8=1TOJ9-1 

1261 A2(J9,J841 )=A2(J9»J8 41 (A20 59 8 )-A2( J9-1 5J8) (4% 8-1) )sNEXTI8 
1262 T1=A2(J9.J9 J-A2( U9 »J9-1 Js IFSGN( T1 )*T1-TL>OTHENNEXTJ9ELSE1264¢ 
1263 PRINT"Tolerance of"3TL3"not met after five extrapolations" 
1264 IN=A2(J9,J9) 

1265 PRINT"Integral ="3:PRINTUSING"S#H###. FHHRH" 5 IN 

1266 PRINT" Actual Tolerance = "3:PRINTUSING'##. ####"5T1*¥SGN(T1) 
1267 SOUND1567 »10:SOUND1244¢ ,10:SOUND1046,10:SOUND783 ,20 

1268 SOUND1046,10:SOUND783 »40 

1269 INPUT"Hit Enter To Continue: "329:GOTO1205 

1275 FORJ7=1T06: FORJ6=1TOJ7:PRINTUSING''##. ###" 5A2( 57,6 )5 

1276 NEXTJ6:PRINT"" :NEXTJ7: INPUTZ9: RETURN 

1280 REM Simpson's Rule Sum 

1281 Y=YU:GOSUB1285:SS=F : Y=YL:GOSUB1285:SS=SS+F 

1282 FORJ5=2TO(N/2): Y=Y+DY :GOSUB1285:SS=SS+4F : Y=Y¥+DY :GOSUB1285 
1283 SS=SS+2*F :NEXTJ5: Y=Y+DY : GOSUB1285: SS =SS+4¢F ;: RETURN 

1285 ‘Flxsy) to be integrated: 

1286 F=1 

1288 'X & Y = independent variables. Hit F8, then FG When Done. 
1289 RETURN 

1290 ‘Limits of Integration: 

1291 ‘XLOWER/XUPPER are constants. 

1292 ‘YUPPER & YLOWER may be constants or given in terms of X. 
1293 XUPPER=1.5707963 

1294 XLOWER=0 

1295 RETURN 

1296 YUPPER=SIN(X ) 

1297 YLOWER=0 

1298 ‘Hit F8, Then F4 When Done 

1299 DY=(YU-YL J/(N+1):RETURN 
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